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Abstract 

Validation  of  the  turbulence  models  in  the  CFD  code  Beggar  for  unsteady  turbulent 
flow  is  discussed.  Six  validation  cases  of  the  code  are  considered,  three  cases  with  the  intent 
of  validating  the  code  without  the  turbulence  model  and  three  cases  to  validate  the  turbu¬ 
lence  model  itself.  The  code  validation  cases  include  a  vortex  in  an  inviscid  atmosphere, 
a  laminar  cylinder,  and  a  laminar  shock-tube.  The  turbulence  model  validation  cases  in¬ 
clude  a  2D  oscillating  airfoil,  a  3D  cylinder  and  a  turbulent  cavity.  Finally,  a  more  realistic 
simulation  of  a  simplified  store  is  examined.  The  turbulence  models  considered  are  the 
Baldwin-Lomax,  Spalart-Allmaras,  and  a  Detached  Eddy  Simulation  (DES)  model.  The 
conclusions  made  deal  with  necessary  prerequisites  to  properly  simulating  unsteady  turbu¬ 
lent  flow  and  model  selection.  The  prerequisites  necessary  in  the  Beggar  code  are  a  second 
order  temporal  discretization  and  the  calculation  of  at  least  three  Newton  dt-iterations  per 
time  step.  The  Baldwin-Lonrax  model  was  found  to  predict  separation  early,  leading  to  a 
breakdown  of  the  flow  structure.  The  Spalart-Allmaras  model  was  found  to  not  properly 
simulate  unsteady  turbulent  flow.  A  recommendation  of  whether  to  use  the  DES  model  was 
not  made  due  to  time  and  computational  constraints  and  temporary  problems  within  the 
Beggar  code. 
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Validation  of  Turbulence  Models 
for  the  Beggar  Code 
in  Unsteady  Flows 

I.  INTRODUCTION 

The  Air  Force  Seek  Eagle  Office  (AFSEO),  Eglin  Air  Force  Base  (AFB),  FL,  is  the 
United  States  Air  Force  (USAF)  authority  for  weapons  certification  efforts.  AFSEO 
performs  test  and  evaluation  for  aircraft/store  compatibility  certification  and  uses  Compu¬ 
tational  Fluid  Dynamics  (CFD)  to  support  this  process.  The  CFD  code  that  the  Computa¬ 
tional  Aeromechanics  Team  (CAT)  has  produced  and  implements  is  Beggar.  Determining 
the  flow  about  an  aircraft/store  combination  can  be  extremely  difficult.  Complicated  ge¬ 
ometry  such  as  pylons,  launchers,  and  internal  weapons  bays  can  create  severe  acoustic  and 
aerothermodynamic  environments,  which  are  challenging  to  simulate  numerically.  The  ad¬ 
ditional  challenge  of  rapidly  and  accurately  simulating  the  trajectory  of  a  store  separation 
in  a  high- volume  simulation  environment  is  beyond  the  capabilities  of  most  CFD  programs. 
Beggar  was  specifically  designed  to  excel  under  these  conditions. 

The  USAF  requirement  for  numerous,  simultaneous  and  quick-reaction  solutions  for 
a  wide  variety  of  stores  and  aircraft  increases  the  need  for  a  fast  and  simple  turbulence 
model.  For  this  reason,  the  CAT  has  used  the  Baldwin-Lomax  (1)  algebraic  model  for  most 
of  its  turbulent  simulations.  Although  the  CAT  believes  that  the  Baldwin-Lomax  model 
has  provided  good  results  for  their  applications,  industry  has  been  moving  towards  more 
sophisticated  models,  such  as  the  Spalart-Allmaras(34)  one-equation  model  and  the  k  —  e 
and  k  —  u)  (46)  two-equation  models. 

After  an  initial  evaluation  of  turbulence  models,  the  CAT  recently  chose  to  include  a 
Spalart-Allmaras  one-equation  model  with  the  option  for  Detached  Eddy  Simulation  (DES) 
in  the  Beggar  code.  This  research  effort  will  address  the  validation  of  the  new  model  and 
will  compare  Spalart-Allmaras  and  Baldwin-Lomax  results  for  a  variety  of  simple  test  cases. 
In  addition,  a  sample  store  will  be  considered  with  both  turbulence  models.  Conclusions  will 
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be  drawn  on  the  abilities  of  the  Baldwin-Lomax  and  Spalart-Allmaras  models  to  properly 
simulate  store  separation  problems. 

As  an  algebraic  model,  Baldwin-Lomax  does  not  include  any  additional  differential 
equations  to  solve  for  the  Reynolds  stresses  (39).  This  makes  it  a  computationally  inex¬ 
pensive  model,  which  is  a  desirable  quality  for  a  production  engineering  code.  A  detailed 
description  can  be  found  in  Section  2.5.  Baldwin-Lomax  has  been  shown  to  properly  simu¬ 
late  flat  plate  boundary  layers(l).  Although  it  has  also  been  noted  that  algebraic  models  are 
unreliable  for  separation  flow(45),  comparisons  with  measurements  from  two  experiments 
show  agreement  in  the  prediction  of  separation  and  reattachnrent  points  within  about  one 
boundary-layer  thickness(l).  Baldwin-Lomax  has  also  proven  itself  excellent  for  attached 
boundary  layers,  however  it  is  not  applicable  to  free-shear  flows. 

Because  it  does  not  include  additional  differential  equations,  the  Baldwin-Lomax 
model  cannot  account  for  turbulent  transport.  In  some  cases,  however,  turbulent  viscosity 
which  exists  upstream  from  a  point  may  significantly  influence  the  solution  at  that  point. 
Without  transport,  the  solution  will  be  poor.  A  simple  example  of  this  is  a  circular  arc 
bump  in  a  two-dimensional  channel.  Baldwin-Lomax  predicts  almost  no  eddy  viscosity  on 
the  bump  where  models  with  transport  show  significant  eddy  viscosity  over  the  bump(19). 
This  problem  may  also  be  present  during  store  separation  when  stores  at  high  angles  of 
attack  experience  flow  separation.  In  some  cases  however,  such  as  attached  boundary  layer 
flows,  Baldwin-Lomax  can  give  as  good  if  not  better  results  than  more  complex  models. 

The  one-equation  Spalart-Allmaras  model  includes  a  transport  equation  derived  using 
empiricism  and  arguments  of  dimensional  analysis  (34).  For  more  details  on  the  implemen¬ 
tation  of  Spalart-Allmaras,  see  reference  (34).  Spalart-Allmaras  is  a  relatively  fast  and 
numerically  forgiving  model.  It  was  made  with  airfoils  in  mind.  For  other  cases,  specifi¬ 
cally  round  jets,  Spalart-Allmaras  may  not  perform  well(33).  The  Spalart-Allmaras  model 
predicts  skin  friction  for  attached  boundary  layers  as  well  as  algebraic  models.  With  the  ex¬ 
ception  of  jets,  Spalart-Allmaras  also  predicts  free-shear  flows  well. (45)  For  attached  flow, 
the  Baldwin-Lomax  model  in  Beggar  should  perform  as  well  as,  if  not  better  than,  the 
Spalart-Allmaras  model.  For  free-shear,  such  as  flow  over  a  cavity,  we  would  expect  the 
Spalart-Allmaras  model  to  outperform  Baldwin-Lomax. 
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Both  the  Baldwin-Lomax  and  Spalart-Allmaras  models  are  Reynolds  Averaged  Navier- 
Stokes  (RANS)  models  which  use  temporal  averaging.  Large  Eddy  Simulation  (LES)  is  an 
alternate  modeling  approach  that  uses  spatial  averaging.  Since  large  eddies  cannot  form  near 
walls,  LES  does  not  perform  well  near  walls  without  refining  the  grid  to  a  Direct  Numerical 
Simulation  (DNS)  scale.  However,  in  the  wake  region,  LES  performs  well.  Baldwin-Lomax 
and  Spalart-Allmaras  perform  well  near  walls,  but  since  they  are  based  on  the  distance 
from  the  wall  they  tend  to  not  perform  well  in  wake  regions.  The  DES  model  attempts 
to  combine  the  best  of  the  two  models.  It  implements  the  RANS  model,  but  adds  spatial 
filtering  to  allow  direct  calculation  of  large  eddies  in  the  wake. 

1.1  Research  Goals 

The  main  goal  of  this  research  is  to  determine  which  turbulence  model,  Baldwin- 
Lomax,  Spalart-Allmaras,  or  DES,  will  perform  best  for  the  CAT.  To  ensure  that  model 
can  function,  the  following  prerequisite  options  will  be  investigated: 

1.  1st  or  2nd  order  temporal  discretization 

2.  Roe  or  Steger- Warming  method  for  calculating  fluxes 

3.  Number  of  Newton  dt-iterations  to  perform  per  time  step 

Additional  considerations  while  evaluating  the  previous  options  are  artificial  dissipation, 
how  well  calculated  data  matches  experiment  or  theoretical  data,  and  cpu-time  required  to 
perform  the  computation.  For  example,  if  the  DES  marginally  outperforms  the  Baldwin- 
Lomax  model,  but  takes  three  times  the  cpu-time,  using  DES  would  not  be  advisable  in  a 
production  code  such  as  Beggar. 

1.2  Prior  Research 

1.2.1  Beggar  History.  Being  a  production  code,  there  has  been  a  fair  amount  al¬ 
ready  written  about  Beggar(4,  6,  12,  11,  14,  26,  28,  29,  30,  42).  Rizk,  Ellison,  and  Prewitt(30) 
give  a  good  overview  of  the  utility  of  Beggar.  Beggar  has  been  designed  primarily  to  cal¬ 
culate  flow  interactions  with  moving  rigid  bodies,  where  the  bodies  may  also  interact.  To 
accomplish  this,  Beggar  uses  an  overset,  or  “Chimera”,  grid  system.  Each  object  has  its 
own  grids  with  its  own  mass.  Some  grid  dynamics  may  be  set  by  the  user.  For  example, 
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a  three  grid  system  includes  a  store  with  a  fin  falling  in  a  large  cartesian  background  grid. 
The  user  may  specify  that  the  fin  grid  will  stay  attached  to  the  store  grid.  Furthermore,  the 
user  may  specify  that  while  staying  attached,  the  fin  grid  will  rotate  about  a  specified  axis, 
to  pitch  the  fin  for  example.  Beggar  then  determines  the  best  points  that  each  interpolation 
cell  should  get  information  from  in  each  overset  grid. 

Westmoreland(42)  discusses  a  comparison  between  using  inviscid  and  viscous  solvers  to 
simulate  store  separations.  Westmoreland  simulated  a  jettison  flight  test  of  the  sample  store. 
He  found  that  when  the  two  solvers  were  run  on  stationary  stores,  substantial  differences  in 
aerodynamic  loads  were  found  in  areas  of  separation  and  boundary  layer  growth.  However, 
when  the  store  was  released,  the  inertial  forces  and  externally  applied  loads  were  two  orders 
of  magnitude  larger  than  the  aerodynamic  forces.  Therefore  the  trajectory  of  the  store  was 
not  influenced  by  the  choice  of  using  an  inviscid  or  viscous  solver.  Noack  and  Jolly(26) 
obtained  similar  results  simulating  the  release  of  a  JDAM  from  an  F-18C  aircraft. 

If  the  major  choice  of  either  an  inviscid  and  viscous  solver  does  not  change  the  tra¬ 
jectory  of  a  store,  there  is  little  chance  that  the  choice  of  which  turbulence  model  to  use 
will  have  any  other  effect  besides  a  dramatic  increase  in  the  required  computing  resources. 
The  case  of  store  trajectory  or  release  is  not  the  only  situation  that  the  AFSEO  simulates, 
however.  There  are  other  cases,  such  as  a  store  falling  from  a  bay,  where  aerodynamic  forces 
will  have  significant  effects,  in  which  this  analysis  have  more  direct  applicability. 

1.2.2  Other  Validation  Efforts.  This  document  closely  follows  the  turbulence 
model  validation  methods  of  Dr  Robert  Nichols(22)  .  Nichols  stresses  that  before  a  turbu¬ 
lence  model  can  be  validated,  the  Navier-Stokes  solver  must  be  tested  to  ensure  that  the 
code  has  the  capability  to  properly  simulate  a  time  accurate  calculation.  Dr  Nichols  places 
large  emphasis  on  ensuring  that  no  solver  component  has  too  much  artificial  dissipation. 
He  has  validated  codes  that  use  a  Newton  iteration  time  step  method  similar  to  Beggar. 
When  validating  these  codes,  many  simulations  have  been  run  while  varying  the  number 
of  Newton  dt-iterations  to  determine  the  minimum  number  required  to  accurately  compute 
a  time  accurate  simulation.  Nichols  concludes  that  to  simulate  a  time  accurate  flow  the 
following  solver  options  are  necessary: 
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1.  Second  order  time  derivatives:  to  lower  numerical  dissipation 

2.  Numerical  flux  algorithms  with  low  dissipation:  to  reduce  numerical  dissipation 

3.  Subiteration  scheme  to  improve  local  convergence:  to  allow  larger  time  steps 

4.  A  rapidly  converging  inner  algorithm:  to  reduce  cpu-time 

Test  cases  to  validate  the  DES  model  have  been  taken  from  references  (23)  and  (21). 

In  reference  (23)  a  two  dimensional  oscillating  pitching  airfoil  and  a  turbulent  cavity  are 
both  examined.  The  oscillating  airfoil  is  a  NACA-0015  oscillating  ±4.2°  around  an  angle  of 
attack  of  a  =  11°  at  a  Reynolds  number  of  Re  =  1.95xl06.  At  these  conditions,  separation 
occurs  during  the  downstroke  cycle.  Using  a  time  step  of  2.2xl0~5  seconds,  Nichols  and 
Nelson  found  that  the  RANS  turbulence  models  did  not  show  any  signs  of  shedding  eddies, 
while  the  DES  model  did  have  shedding  eddies.  Similar  effects  were  found  in  the  turbulent 
cavity,  where  much  more  unsteady  motion  was  apparent  with  the  DES  model  than  with  the 
RANS  model. 

The  test  case  taken  from  reference  (21)  is  a  3D  circular  cylinder  in  cross-flow.  Nichols 
attempts  a  grid  convergence  study,  but  notes  that  the  finer  a  grid  is,  the  smaller  the  turbulent 
length  scales  are  which  can  be  captured  by  the  grid.  This  continues  until  the  smallest  of 
the  turbulent  length  scales  is  reached.  Currently,  it  is  not  feasible  computationally  to  arrive 
at  grids  this  fine.  Therefore,  a  true  grid  convergence  study  could  not  have  been  performed. 
To  attempt  to  check  grid  convergence,  the  average  drag  coefficient  was  recorded  for  three 
cylinder  grids  of  varying  levels  of  refinement.  Nichols  found  that  the  DES  model  did  not 
reach  grid  convergence  with  the  three  grids.  In  his  temporal  study,  Nichols  found  that  200 
time  steps  per  shedding  frequency  are  necessary  for  temporal  accuracy  using  second  order 
temporal  discretization  and  a  hybrid  RANS/LES  model.  Similar  to  the  turbulent  cavity 
and  oscillating  airfoil,  Nichols  found  that  the  LES  model  showed  weak  eddies  in  the  far  wake 
where  the  RANS  models  damped  them  out. 

In  another  study,  Nichols  and  Nelson(18)  studied  the  effectiveness  of  a  hybrid  RANS/LES 
model  they  developed  and  made  comparisons  to  the  DES  model.  Their  model  implemented 
a  k  —  e  model  to  account  for  subgrid  scale  turbulence.  The  test  case  in  (18)  was  a  compress¬ 
ible  mixing  layer  where  one  layer  was  subsonic  and  the  other  was  supersonic.  Nichols  and 
Nelson  found  that  both  the  Nichols-Nelson  and  DES  models  showed  drastic  improvements 
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over  the  RANS  models.  Their  results  matched  well  with  traditional  LES  models.  They 
would  not  make  further  conclusions,  as  no  single  source  of  error  could  be  identified. 

Seror,  Rubin,  Peigin,  and  Epstein(32)  discussed  the  validation  of  a  Spalart-Allmaras 
turbulence  model  in  a  CFD  code  called  “NES”.  Their  paper  provides  insight  into  validat¬ 
ing  the  Spalart-Allmaras  model  without  being  overshadowed  by  the  validation  of  the  DES 
model.  Similar  to  Beggar,  the  default  turbulence  model  in  NES  was  Baldwin-Lomax.  Due 
to  the  difficulty  in  the  fine  grid  resolution  required  by  Baldwin-Lomax  and  the  poor  results 
it  gave,  Seror  et  al.  decided  to  include  the  Spalart-Allmaras  model.  The  test  case  they 
ran  was  a  2D  RAE2822  supercritical  airfoil  at  transonic  flight  conditions.  Seror  et  al  found 
that  the  Spalart-Allmaras  model  more  accurately  predicted  the  shock  location  than  the 
Baldwin-Lomax  model,  and  therefore  improved  the  prediction  of  lift  and  drag.  They  also 
noted  that  there  was  much  variation  in  the  wake  between  the  two  models. 


1.3  Research  Approach 

Validation  of  the  additional  turbulence  models  is  broken  down  into  two  main  com¬ 
ponents:  a  non-turbulent  study  and  a  turbulent  study.  The  individual  cases  are  shown  in 
Table  1. 


Table  1:  Test  cases 


Case 

Study 

Convecting  Vortex 

Non- Turbulent 

Shock-Tube 

Non- Turbulent 

2D  Vortex  Shedding  Cylinder 

Non- Turbulent 

2D  Oscillating  Airfoil 

Turbulent 

3D  Vortex  Shedding  Cylinder 

Turbulent 

Turbulent  Cavity 

Turbulent 

In  addition  to  these  problems,  the  turbulence  models  will  be  applied  to  a  store  sep¬ 
aration  problem.  Results  will  be  used  to  determine  the  effects  of  using  one-equation  or 
algebraic  turbulence  models  on  simulating  store  separation. 

1.3.1  Non- Turbulent  Study.  The  non-turbulent  study  includes  three  test  cases 
including  a  convecting  vortex,  a  shock-tube,  and  a  2D  vortex  shedding  cylinder.  These 
cases  are  presented  in  order  of  difficulty. 
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The  first  case  examined  in  the  process  of  validating  the  turbulence  model  is  a  vortex 
convecting  in  an  inviscid  atmosphere.  The  inviscid  convecting  vortex  is  a  quick  simulation 
to  run  to  check  the  numerical  dissipation  in  each  component  of  the  code. 

The  shock-tube  and  2D  cylinder  cases  are  simulated  with  the  intention  of  determining 
the  proper  number  of  Newton  dt-iterations  to  use.  Both  of  these  cases  are  quick  to  solve, 
allowing  many  different  solver  option  combinations  to  be  simulated. 

1.3.2  Turbulent  Study.  The  turbulent  study  has  three  test  cases  consisting  of  a 
2D  oscillating  airfoil,  a  3D  vortex  shedding  cylinder,  and  a  turbulent  cavity. 

The  first  turbulent  case  is  a  simple  oscillating  2D  airfoil.  Two  initial  angle’s  of  attack, 
a0,  are  examined,  a0  =  4.0°  and  aQ  =  11°.  With  the  aQ  =  4°,  the  angle  of  attack  is  always 
small  enough  that  the  flow  remains  attached  to  the  airfoil  throughout  the  simulation.  This 
case  will  be  used  as  another  check  to  see  that  the  solver  can  accurately  model  a  time  accurate 
simulation.  The  larger  aD  will  have  separation  off  the  top  of  the  airfoil,  which  will  provide 
insight  into  the  turbulence  models  effect  on  separated  flow. 

A  vortex  shedding  turbulent  cylinder  is  more  difficult  than  the  2D  case.  With  Re  > 
180,  the  flow  can  no  longer  be  approximated  as  2D.  Numerical  data  will  be  compared  to 
experimental  data. 

An  unsteady  simulation  of  the  Weapons  Internal  Carriage  and  Separation  (WICS) 
weapons  bay  will  also  be  simulated.  Experimental  data  has  been  gathered  at  two  specific 
points  inside  the  cavity.  Care  will  be  taken  to  ensure  that  the  CFD  grid  will  have  a  data 
point  that  matches  where  the  experimental  data  was  taken.  CFD  data  will  then  be  compared 
to  the  experimental  data.  A  grid  convergence  study  will  be  performed  for  this  case. 

1.3.3  Store  Separation  Study.  Although  important  for  validation,  the  previous 
cases  are  all  very  simple,  and  do  not  apply  directly  to  everyday  production  simulations. 
To  have  some  idea  of  the  effects  that  the  different  turbulence  models  will  have  on  a  more 
realistic  case,  a  sample  store  jettison  will  be  simulated.  As  the  turbulence  model  will  have 
the  most  effect  on  the  store  when  inertial  effects  are  minimal,  the  initial  motion  of  a  store 
falling  from  rest  will  be  considered  separate  from  when  the  store  has  significant  inertia.  The 
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time  step  will  then  be  gradually  increased  to  allow  the  store  to  move  a  noticeable  distance 
in  a  reasonable  amount  of  time  steps. 

1.4  Document  Organization 

First  a  background  will  be  presented  to  familiarize  the  reader  with  the  terminology 
used  in  this  document.  Chapter  II  will  start  with  the  governing  equations  of  aerodynamic 
flow.  Then  a  description  of  how  Beggar  implements  these  equations  is  given.  A  brief 
introduction  to  turbulence  follows,  ending  with  how  turbulence  is  modeled  in  CFD  codes. 
Chapter  III  is  dedicated  to  checking  Beggar  to  ensure  it  has  the  capability  of  properly 
simulating  time  accurate  flow.  Results  will  be  given  along  with  the  presentation  of  the  case. 
When  that  has  been  done,  the  turbulence  validation  cases  will  be  presented  in  Chapter  IV. 
Again,  results  will  be  included  with  the  presentation  of  the  case.  Chapter  V  will  present  the 
sample  store  jettison  simulation  and  results.  Finally,  some  conclusions  and  recommendations 
will  be  given  in  Chapter  VI. 


II.  BACKGROUND  AND  THEORY 


2.1  General  Aerodynamic  Flow 

For  most  of  the  test  cases  used  to  validate  the  new  turbulence  models  in  the  Beggar 
code,  the  Navier-Stokes  equations  will  be  used  to  solve  the  flow  physics. 

Historically  speaking,  the  Navier-Stokes  equations  are  the  viscous  conservation  of  mo¬ 
mentum  equations  derived  from  Newton’s  second  law  and  the  deformation  law  for  a  new- 
tonian  fluid (43).  However,  in  CFD  the  name  Navier-Stokes  is  given  to  the  set  of  equations 
including  the  conservation  of  mass,  momentum  and  energy.  The  Beggar  code  is  a  finite 
volume  solver,  which  uses  the  integral  form  of  the  Navier-Stokes  equations: 


J  ^dV  +  J(FC  -  Fv)dA  =  0  (1) 
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is  the  work  done  by  the  viscous  stresses  and  the  heat  conduction  of  the  fluid(2),  U  =  V  •  h 
is  the  contravariant  velocity,  and  r  is  the  stress  tensor  defined  from  the  deformation  law  for 
a  Newtonian  fluid  given  by: 
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which,  using  Stokes  hypothesis  (A  =  —  |/i),  reduces  to: 
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Equation  4  is  valid  for  fluids,  including  most  gases,  with  the  following  properties  as  noted 
by  Stokes (43): 


1.  There  is  a  linear  relationship  between  stress  and  strain. 

2.  An  isotropic  fluid  is  assumed. 

3.  As  strain  rates  approach  zero,  the  shear  stresses  will  approach  zero  and  the  normal 
stresses  will  equal  the  static  pressure 

4.  Body  forces  are  absent. 

While  Beggar  uses  dimensional  quantities  to  calculate  body-to-body  and  flow-to-body 
interactions,  the  flow  solver  is  non-dimensional.  To  non-dimensionalize  the  Navier-Stokes 
equations,  the  values  shown  in  Equation  5  are  inserted  into  Equation  1.  The  non-dimensional 
quantities  in  Equation  5  are  represented  with  an  asterisk. 


P  —  P  Poo 


Et  =  E\  p 00*^00  P  =  P*  Poodle 


t  =  t 


:  ETef 


(5) 


u  =  u*aoo  v  =  v*aco  iu  =  w*a00 


Situations  exist  where  viscosity  does  not  play  a  significant  role  in  determining  the 
flow  properties.  For  example,  in  high  Reynolds  number  flows,  the  inertial  effects  outweigh 
the  viscous  effects.  In  these  cases  viscosity  is  only  important  in  a  thin  layer  that  lies  on 
body  surfaces.  The  region  outside  of  this  boundary  layer  may  be  treated  as  inviscid  with  no 
spatial  temperature  gradient  effects.  For  these  situations,  the  Navier-Stokes  equations  may 
be  simplified  by  eliminating  all  viscous  terms,  which  results  in  the  Euler  equations.  The 
Euler  equations  are: 


IMdv+l/<iA=o  <6) 

where  Fc  and  Q  are  the  same  as  in  Equation  1.  The  convecting  vortex  case  is  solved  with  the 
Euler  equations  in  place  of  the  Navier-Stokes  equations.  The  fact  that  the  Euler  equations 
do  not  include  any  effects  from  viscosity  will  be  exploited  to  help  determine  which  CFD 
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solver  options  to  use.  Further  discussion  of  the  convecting  vortex  case  can  be  found  in 
Section  3.1. 

The  equations  described  above  are  the  basis  of  general  aerodynamic  flow.  The  Navier- 
Stokes  equations  include  everything  needed  to  solve  turbulent  flow.  However,  to  solve  for 
turbulent  flow  using  only  the  Navier-Stokes  equations  is  very  problematic  and  is  rarely  done. 
This  is  described  in  Section  2.4. 


2.2  Beggar  Implementation 

The  Navier-Stokes  equations  (Equation  1)  are  continuous  and  can  be  solved  at  all 
points  in  a  flow  at  all  times.  However,  an  exact  general  solution  to  the  Navier-Stokes 
equations  has  not  been  found  and  is  not  expected  to  be  found  any  time  soon.  Therefore, 
they  are  solved  numerically.  To  do  this,  the  Navier-Stokes  equations  are  discretized  to  be 
solved  at  specific  times  and  locations.  An  implicit  discretization  of  Equation  1  is  given  by: 


-v  +  ^(J^+1-r„"+1)|AI=o 


ann+1 

where  —  is  the  time  discretization  of  the  the  change  in  the  conserved  variables  Q  in  the 
cell  volume  V  with  respect  to  time  and  the  second  term  is  the  sum  of  the  fluxes  through  the 
sides  of  the  cell.  Superscript  n  is  the  specific  time  step.  It  is  difficult  to  compute  the  fluxes 
at  time  n+  1.  To  simplify  the  computation,  the  fluxes  are  linearized  in  time  which  results 
in(15): 

F?+1  =  F?  + -^AQn+1  (8) 


Fy+]  =  Fy  +  ~0^AQn+] 

Substituting  equations  8  and  9  into  equation  7  gives 


dQn+1  v\^(dFc  dFv 
dt  +^\dQ  dQ 


)  AQn+1\A\  +  E  (Ecn  -  F?)  \A\  =  0 
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Pulling  |  A  |  into  F  and  simplifying  gives 


dU 

dQ 


dQn+1 

~dt~ 


A  Qn+1  =  -Kn 


(11) 


where  is  the  time  discretization  of  the  conserved  variables,  is  the  implicit  oper¬ 

ator,  or  Jacobian,  which  is  derived  using  first  order  Steger- Warming  fluxes  (30),  and  1Zn  is 
the  explicit  term  commonly  referred  to  as  the  “right  hand  side”  (RHS) 

Two  different  time  discretizations  are  compared  in  this  study,  the  first  order  backward 
Euler  time  discretization  and  the  second  order  three  point  backward  time  discretization.  The 
backward  Euler  time  discretization  method  is  given  as(30): 


dQn+1  _  Qn+1  -  Qn 
dt  At 


(12) 


and  the  three  point  backward  time(8)  discretization  is  given  as: 


dQn+1  3 Qn+1  -  4 Qn  +  Q 


n—  1 


dt 


2A  t 


(13) 


The  three  point  backwards  time  discretization  is  a  higher  order  and  generally  more  accurate 
method  than  the  backward  Euler  time,  but  it  requires  storing  all  of  the  conserved  variables 
for  the  two  previous  time  steps  whereas  the  backward  Euler  time  method  only  requires 
storing  the  previous  time  step. 

Two  different  methods  for  calculating  the  right  hand  side  are  included  in  Beggar, 
and  each  will  be  tested;  the  second  order  Roe(31)  and  second  order  Steger- Warming(36) 
methods.  Both  of  these  methods  are  flux  splitting  methods.  The  Steger- Warming  method 
splits  the  fluxes  into  negative  and  positive  contributions  based  on  eigenvalues.  The  Roe 
method (2)  is  an  approximate  Riemann  solver  that  splits  the  difference  of  the  fluxes. 

Accurately  calculating  each  time  step  is  critical  to  unsteady  flows.  To  ensure  that  each 
time  step  has  been  computed  accurately,  Beggar  solves  a  user  specified  number  of  Newton 
dt-iterations  for  each  time  step.  Newton’s  method(44)  for  a  vector  function  G(x)  =  0  is 
given  as: 

G'{xm)(xm+1  -  xm )  =  —G(xm)  (14) 
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where  superscript  m  is  the  Newton  dt-iteration  and  G'(x)  is  the  Jacobian  matrix  given  by: 


an  (a)  ai2(x)  ...  a1;fe(x) 

a2i(x)  a2  2{x)  ...  a2,k(x) 


G\x)  = 


Qkii.x')  o,k2{x)  ^k,k(.x) 


where 


a>ij(x)  = 


dGj{x) 

dxj 


Rearranging  Equation  11  to  resemble  Equation  14  leads  to: 


fdR 

\dQ 


n,m 


A  Q 


n+1, m+1 


'  dQn+1’m 

di 


+  Rn’m 


(15) 


(16) 


(17) 


where 


A  Q 


n+1, m+1 


Qn+l,m+l 


Qn+l,m 


(18) 


To  solve  for  Qn+1,m+1 ,  the  AQn+1,m+1  term  must  be  isolated.  This  would  involve  inverting 
the  Jacobian,  whose  dimensions  may  be  on  the  order  of  millions.  Inverting  a  matrix  this 
large  is  intractible.  Beggar  avoids  inverting  the  Jacobian  by  solving  for  Qn+1>m+1  iteratively 
using  the  Symmetric  Gauss-Seidel  method.  The  Gauss-Seidel  method  solves  the  equation: 


[A]x  =  b 


(19) 


for  x  iteratively  by  breaking  [A]  into  an  upper  and  a  lower  triangle  and  solving 


x\+1  =  bi-  I+x^1)  -  C+xfc) 


fc+B 


(20) 


where 


L*(xfc+1)  =  y~]  Aitjx)+1 
j= i 


(21) 
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and 


Jmax 

Ui{x.k)  =  (22) 

j=i+ 1 

where  k  +  1  implies  an  updated  value.  Alternating  starting  from  the  first  element  and  the 
last  element  makes  the  process  “symmetric.”  The  Symmetric  Gauss-Seidel  iterations  are 
hereafter  referred  to  as  inner  iterations.  Equation  20  is  solved  iteratively  until  convergence 
criteria  is  met  or  the  user  set  maximum  number  of  Symmetric  Gauss-Seidel  iterations  is 
computed. 

2. 3  Turbulence 

Before  turbulence  modeling  is  discussed,  a  brief  introduction  to  turbulence  is  pre¬ 
sented.  Turbulence  is  an  important  aspect  of  computational  fluid  dynamics  since  most 
engineering  flows  are  turbulent.  In  non-turbulent  flows,  random  fluctuations  in  the  flow  are 
damped  out  by  viscous  forces.  In  turbulent  flow,  the  inertial  forces  overpower  the  viscous 
forces  and  the  random  fluctuations  are  no  longer  damped  out,  and  actually  grow.  These 
fluctuations  become  an  important  characteristic  of  the  flow,  communicating  with  neighbor¬ 
ing  fluctuations  and  creating  vortices  of  different  scales,  which  also  communicate  with  each 
other  via  an  energy  cascade. 

Figure  1  shows  a  typical  time  history  of  velocity  measured  at  a  single  point  in  a 
turbulent  flow.  It  illustrates  the  different  scales  that  are  present  in  a  turbulent  flow.  There 
is  a  lot  of  high  frequency  turbulence  which  appears  as  noise  in  the  figure.  If  the  sampling  rate 
was  reduced  to  filter  out  the  high  frequency  turbulence,  the  general  shape  of  the  curve  would 
still  remain.  This  shape  would  be  the  largest  turbulent  scale,  known  as  the  integral  scale. 
By  sampling  more  frequently,  smaller  patterns  will  appear  on  the  larger  shape,  similar  to 
zooming  in  on  a  fractal.  There  is  a  limit  when  sampling  more  frequently  will  not  produce  any 
more  patterns  and  the  curve  will  smooth  out.  At  this  sampling  rate  the  smallest  turbulent 
scale,  known  as  the  Kolmogorov  scale,  can  be  determined. 

Figure  2  displays  important  characteristics  of  the  different  length  scales.  On  the 
abscissa  is  the  frequency,  /,  which  is  inversely  proportional  to  wavelength.  Therefore,  the 
smaller  the  abscissa  value,  the  larger  the  length  scale.  From  largest  to  smallest,  the  three 
main  turbulent  length  scales  are  the  integral  scale,  the  Taylor  scale,  and  the  Kolmogorov 
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Velocity 


Figure  1:  Velocity  with  respect  to  time  at  a  point  in  turbulent  flow. 


PSD 


Figure  2:  Turbulent  energy  spectrum. 
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scale.  Each  scale  appears  in  its  own  range  in  Figure  2.  The  integral  scale  is  where  turbulence 
obtains  its  energy.  The  integral  scale  interacts  with  the  mean  shear  and  gains  kinetic  energy 
from  it. 

Following  the  integral  scale  is  the  Taylor  scale,  which  resides  in  the  middle  range 
known  as  the  inertial  range.  In  this  range  energy  is  transferred  from  larger  scale  eddies 
to  smaller  scale  eddies  without  significant  losses  from  viscosity.  Using  known  information 
about  the  small  scale  and  the  large  scale  and  performing  a  dimensional  analysis, (16)  Andrei 
Kolmogorov  derived  the  Kolmogorov  power  law, 

E(X,t)  =  Cel /“I  (23) 

where  E  is  energy,  e  is  dissipation  and  /  is  frequency.  As  seen  in  Figure  2,  the  inertial 
range  has  a  constant  slope  with  a  value  of  /_s  on  an  energy  spectrum  plot,  verifying  the 
Kolmogorov  power  law.  The  Kolmogorov  power  law  has  been  verified  so  well  that  the 
inertial  range  is  defined  as  the  region  on  a  power  spectrum  density  plot  with  a  slope  of  —  | . 

The  third  scale  has  been  named  after  Kolmogorov.  The  Kolmogorov  scale  is  the 
smallest  scale.  At  this  size,  the  eddies  are  too  small  to  have  enough  inertia  to  overcome 
viscosity.  Energy  is  then  transferred  from  kinetic  energy  to  heat  via  viscous  dissipation. 
The  Kolmogorov  scale  is  dependent  on  two  things,  the  energy  it  receives  from  the  larger 
eddies  and  the  viscosity  of  the  fluid.  Combining  these  two  variables  in  such  a  way  that 
results  in  units  of  length  gives  the  Kolmogorov  length  scale,  rj,  where 

V  =  ( —  J  oc  Ret  4  (24) 

and 

Ret  =  v!  L/v  (25) 

The  Kolmogorov  scale  appears  at  the  far  right  of  the  power  spectrum  density  plot  shown  in 
Figure  2. 

The  Navier-Stokes  equations  contain  everything  needed  to  solve  over  all  ranges  and 
scales  of  a  turbulent  flow.  The  process  of  solving  turbulence  using  nothing  but  the  Navier- 
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Stokes  equations  is  referred  to  as  Direct  Numerical  Simulation  (DNS).  A  DNS  solution 
requires  a  mesh  that  is  fine  enough  to  capture  the  smallest  length  scales  and  a  time  step  small 
enough  to  capture  its  frequency.  Equation  24  shows  that  the  Kolmogorov  scale  gets  smaller 
as  the  Reynolds  number  grows.  Even  at  moderate  Reynolds  numbers  the  Kolmogorov  scale 
is  so  small  that  only  the  simplest  of  cases  can  be  simulated.  For  this  reason,  turbulence 
modeling  is  used. 

2.4  Turbulence  Modeling 

Turbulence  is  a  chaotic  and  complex  phenomenon.  Fluctuations  occurring  in  the  flow 
do  not  follow  a  set  pattern.  For  this  reason,  turbulence  is  best  described  mathematically 
with  statistics.  One  approach  is  to  model  the  flow  in  terms  of  the  sum  of  a  time  averaged 
component  and  a  fluctuating  component,  as  shown  in  Equation  26.  An  overbar  denotes  a 
time  averaged  value  while  a  fluctuation  component  is  denoted  with  an  apostrophe. 

u  =  u  +  vl  p  =  p  +  p' 

v  =  v  +  v'  p  =  p  +  p'  (26) 

w  =  w  +  w'  T  =  T  +  T’ 

Equation  27  shows  rules  for  averaging.  Inserting  the  values  in  Equation  26  into  Equa¬ 
tion  1  and  averaging  the  resulting  expression  using  the  rules  in  Equation  27  results  in  the 
RANS  equations.  The  RANS  equations  assume  that  turbulent  time  scales  are  much  smaller 
than  the  characteristic  time  of  the  flow.  Since  this  study  will  be  using  RANS  models  to 
simulate  unsteady  flow,  care  must  be  taken  to  ensure  that  the  unsteady  characteristic  time 
is  orders  of  magnitude  larger  than  the  averaging  time.  Otherwise,  the  unsteady  effects  will 
be  averaged  out. 


/'  =  0  /  =  /  fg  =  fg  fg  =  0 

_  _  (27) 

f+g=l+g  %  =  %  fg  =  Jg  +  f'g'  ffds  =  fjds 

The  RANS  equations  can  be  simplified  by  applying  Morkovin’s  hypothesis(2),  which 
states  that  the  fluctuating  components  of  density  are  much  smaller  than  the  average  density, 
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p'  «  p.  When  Morkovin’s  hypothesis  is  applied  to  the  RANS  equations  all  the  fluctuating 
terms  except  one  drop  out  of  the  equations.  This  new  term,  shown  in  Equation  28,  is 
called  the  Reynolds  stress  tensor.  The  Reynolds  stress  tensor  captures  the  influence  of  the 
turbulent  fluctuations  on  the  mean  flow.  Being  a  symmetric  tensor,  the  Reynolds  stress 
tensor  adds  6  new  unknowns  to  the  Navier-Stokes  equations. 


^  pu'u' 

pu'v' 

pu'w'  ^ 

(28) 

pu'v' 

pv'v' 

pv'w' 

y  pu'w' 

pv'w' 

pw'w'  j 

A  further  simplification  of  the  RANS  equations  is  the  Boussinesq  assumption(l).  Boussinesq 
assumed  that  the  Reynolds  stress  is  similar  to  the  laminar  viscous  stresses  and  may  be 
written  as 

-pu'jUj  =  2  mSij  (29) 


where 

o.  =  1  + 

l}  2  \dxj  dxi ) 

is  the  strain  rate  tensor.  Therefore,  the  Reynolds  stress  tensor  is  reduced  to  computing  one 
new  variable,  the  turbulent  viscosity  pt-  Now,  the  same  Navier-Stokes  equations  as  used  for 
laminar  flow  shown  in  Equation  1  may  be  used  for  the  turbulent  simulations  with  the  RANS 
turbulence  models,  except  the  p  in  Equation  4  is  replaced  with  the  sum  of  the  laminar  and 
turbulent  viscosities,  p  =  pi  +  pt-  Finding  a  value  for  pt  is  the  goal  of  the  turbulence 
model.  The  two  turbulence  models  under  consideration  in  Beggar,  the  Spalart-Allmaras 
and  Baldwin-Lomax  models,  both  use  the  Boussinesq  assumption. 


2.5  Baldwin-Lomax 

The  Baldwin-Lomax  model(l)  is  a  “zero-equation”  model  which  provides  an  algebraic 
closure  for  the  Reynolds  Averaged  Navier-Stokes  (RANS)  equations  based  on  the  Boussinesq 
assumption.  Baldwin-Lomax  uses  Prandtl’s(43)  mixing  length  theory,  which  by  analogy 
with  kinetic  theory  relates  the  turbulent  viscosity  to  some  turbulent  length  scale  and  a 
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velocity  gradient: 


(31) 


,9  du 

m  05  P‘si 

where  l  is  the  “mixing- length” .  Since  turbulence  is  a  continuum  phenomena(43) ,  the  mixing- 
length  is  much  larger  than  the  mean  free  path  of  the  fluid  particle. 


The  Baldwin-Lomax  model  breaks  the  spectrum  of  turbulent  eddies  into  two  regions, 
the  inner  “near  wall”  region  and  an  outer  “wake”  region.  The  viscosity  in  the  inner  region 
is  calculated  with  a  formula  based  on  Prandtl’s  mixing  length  theory  known  as  the  Prandtl- 
Van  Driest(l)  formula  for  the  inner  layer, 


(dinner  =  P?  Iw 


(32) 


where  to  is  the  magnitude  of  vorticity 


I  /  du  dv\ 2  f  dv  dw\ 2  f  dw  du\2 

U  y  \ dy  dx )  \ dz  dy  )  +  \dx  dz) 

Prandtl  and  Karman(43)  found  that  within  the  inner  layer,  the  mixing-length  was  strongly 
dependent  on  the  distance  to  the  wall  according  to  the  following 


y 2  in  the  viscous  sub-layer  ( y+  <  5) 


l  oc  < 


ky  in  the  overlap  layer  ( y+  <  300) 


(34) 


where 


_i_  v  PwTw 

y  = - y 

pw 


(35) 


Using  this  knowledge,  a  curve  was  generated  that  had  the  form  of  l  oc  y  with  a  damping 
function  added  to  fit  the  parabolic  curve  for  y+  <  5.  This  inner  mixing-length  is  given  by 


l  =  ky  ^1  —  e  y+/A+^j 


(36) 


where  k  is  the  Karrnan  constant  given  as(9) 


k  =  0.41 
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and 


A+  =  26 


which  is  optimized  for  zero  pressure  gradient  simulations. 

The  turbulent  viscosity  in  the  outer  region  is  given  by 

(^i)  outer  F  CcppFwdke^Klebiy)  (37) 


where  the  Clauser  constant  is 


K  =  0.0168 


and 

Ccp  =  1.6 


Fxuake  —  min 


(y 


max  Fmax  i  F-w  k  l-jrnax  ^  dif  /  Fmax  ) 


(38) 


Cwk  =  1.0 


Baldwin  and  Lomax  used  a  value  of  Cwk  =  0.25,  but  Swanson  and  Turkel  (38)  found  that 
this  value  performs  poorly  in  transonic  flows  over  an  airfoil.  Swanson  and  Turkel  recommend 
using  Cwk  =  1.0.  Fmax  and  ymax  are  determined  from 


F(y)  =y\u\(l-e  y+/A+^j 


(39) 


with  Fmax  being  the  largest  value  in  a  particular  profile  (constant  stream-wise  location) 
and  ymax  being  the  cross-flow  location  of  that  maximum.  In  wakes  the  exponential  term  in 
Equation  39  is  set  equal  to  zero(l).  The  Klebanoff  intermittence  factor  is 


FlCleb{y ) 


1  +  5.5 


F KlebV 
Vinax 


(40) 


This  provides  a  correction  for  the  intermittence  in  the  boundary  between  the  turbulent 
and  laminar  flow  regions.  Baldwin  and  Lomax  optimized  Cxieb  and  agreed  that  C-Kieb  =  0.3 
was  the  best  value  for  separated  flow(l).  It  is  common  practice  to  use  Ccp  =  1.6.  The 
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minimum  velocity  magnitude: 


Udif  ~  I  max 


-IVI 


(41) 


is  zero  except  in  wakes. 

Baldwin-Lomax  is  a  simple  to  implement  model  that  provides  good  results  in  attached 
turbulent  flow(45).  However,  there  is  no  turbulent  transport,  meaning  that  the  upwind 
turbulence  is  not  felt  downwind.  A  turbulence  model  which  includes  turbulent  transport  is 
the  Spalart-Allmaras  one-equation  model. 

2.6  Spalart-Allmaras 

The  Spalart-Allmaras  turbulence  model  is  a  one-equation  model,  meaning  that  it 
adds  one  partial  differential  equation  that  must  be  solved  along  with  the  Navier-Stokes 
equations.  Along  with  the  new  equation  comes  another  variable  that  must  be  computed. 
The  new  variable  in  the  Spalart-Allmaras  model  is  v,  a  working  variable  with  units  equivalent 
to  that  of  viscosity.  The  new  equation  introduced  by  Spalart  and  Allmaras,  like  most 
transport  models,  is  similar  to  the  species  equation  for  reacting  flow(34).  The  Spalart- 
Allmaras  model  was  derived  from  empiricism,  dimensional  analysis,  and  Galilean  invariance. 
It  was  developed  with  wall  bounded  flows  in  mind(7).  The  specific  equation  solved  in  the 
Spalart-Allmaras  model  is: 

%  +  Ui^.  =  \  [V  •  (("  +  +  Cfe2(W)2]  +  P{y)  -  D{y)  (42) 

After  solving  Equation  42,  the  turbulent  viscosity  pt  may  then  be  found  from 

IH  =  pvfviix)  (43) 

where 

X  =  —  and  p  =  pu 
v 
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The  viscous  damping  function,  fv\,  is  given  by: 


fv  1  = 


X 


X3  +  C®!  ’ 

P[y)  in  Equation  42  is  the  production  term,  given  by: 


P{v)  -  Cb  1  (  S+  -^pfv'2 


where  S  is  the  magnitude  of  vorticity  and 


(44) 


(45) 


The  destruction  term  in  Equation  42  ,  D(y ),  is  given  as 

D(D)  =  Cwlfw  (^j 


where 


and 


fw  —  9 


1  +  ct 


w  3 


C6W  3 


1/6 


g  =  r+  Cw2(re  -  r) 


(46) 


(47) 


(48) 


(49) 


and 

(Sn2d2  +  vfv2) 

In  both  Equation  45  and  Equation  47,  d  is  the  distance  to  the  closest  wall.  The  remaining 
unknowns  in  the  Spalart-Allmaras  model  are  constants  which  were  determined  numerically 
or  via  matching  empirical  data,  as  shown  below. 


cbl 

=  0.1355 

Cb2  =  0.622 

2 

G  —  — 

3 

Cvl  =  7.1 

Cw  1 

_  Cbi  (i  +  Cb2) 
k2  a 

Cw  2  =  0.3 

Cw  3  =  2.0 

k  =  0.41 
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2.7  Detached  Eddy  Simulation 

The  third  and  final  turbulence  model  considered  in  this  document  is  the  Detached 
Eddy  Simulation  (DES)  model.  DES  is  a  hybrid  of  the  Spalart-Allmaras  model  combined 
with  similar  spatial  filtering  styles  of  Large  Eddy  Simulation  (LES).  All  temporally  aver¬ 
aged  RANS  turbulence  models,  include  Baldwin-Lonrax  and  Spalart-Allmaras,  place  the 
responsibility  of  capturing  turbulence  completely  on  the  model.  All  scales  of  turbulence, 
from  the  smallest  to  the  largest,  must  be  accounted  for  by  these  models.  LES  models,  on  the 
other  hand,  break  the  turbulent  eddies  into  two  scales,  a  grid  realized  region  and  a  subgrid 
region.  The  larger  scales  of  turbulence,  which  are  able  to  exist  away  from  solid  boundaries, 
are  usually  larger  than  the  grid  spacing.  As  mentioned  before,  the  Navier-Stokes  equa¬ 
tions  include  everything  necessary  to  simulate  turbulence,  as  long  as  the  grid  spacing  is  fine 
enough  to  capture  the  eddies  and  the  time  step  is  small  enough  to  capture  the  turbulent 
frequencies.  The  large  turbulent  eddies  away  from  a  wall  are  often  large  enough  that  they 
will  be  captured  by  a  reasonably  fine  grid.  Therefore,  applying  a  RANS  turbulence  model 
will  compute  the  energy  of  the  grid  resolved  eddies  twice,  once  from  the  Navier-Stokes  equa¬ 
tions  and  again  with  the  turbulence  model.  To  correct  for  this,  the  DES  model  combines 
the  Spalart-Allmaras  RANS  model  with  a  spatial  filtering  function.  As  seen  in  Equation  47, 
the  turbulent  destruction  term  of  the  Spalart-Allmaras  model  is  inversely  proportional  to 
the  square  of  the  distance  to  the  nearest  wall,  d.  Away  from  a  wall,  this  term  will  approach 
zero  while  the  production  term  will  approach  a  constant  non-zero  value.  To  account  for 
turbulent  destruction  away  from  a  wall,  the  DES  model  replaces  d  in  the  Spalart-Allmaras 
model  equations  with  d  where 

d  =  min(d,  Cdes  A)  (51) 

where  d  is  still  the  distance  to  the  nearest  wall,  A  is  the  largest  grid  spacing  of  the  cell,  and 
Cdes  is  an  empirically  determined  constant  set  to  Cdes  =  0.65. 

2.7.1  Beggar  Implementation.  The  turbulence  models  in  Beggar  are  solved  along 
with  the  flow  equations.  Therefore,  the  turbulent  variable  is  updated  every  Newton  dt- 
iteration.  Then  the  flow  equations  are  solved.  Baldwin-Lonrax  is  an  algebraic  model,  so  it 
explicitly  calculates  and  updates  the  turbulent  variable  each  time  step.  Because  Spalart- 
Allmaras  solves  a  differential  equation,  it  is  broken  up  into  Jacobian  and  right  hand  side 


23 


term  like  the  flow  equations  and  solved  in  the  same  manner  as  the  flow  equations.  The  DES 
model  is  solved  in  the  same  manner  as  Spalart-Allmaras. 
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III.  NON-TURBULENT  STUDY 


There  are  certain  prerequisites  that  any  CFD  simulation  must  meet  before  any  turbulence 
model  can  perform  well.  The  code  cannot  have  overwhelming  artificial  dissipation,  or  the 
effects  of  turbulent  viscosity  will  be  masked.  Time  steps  and  grid  cells  must  not  be  too 
large,  or  the  steep  gradient  will  not  be  captured.  Finally,  each  time  step  must  be  sufficiently 
converged,  or  time  accuracy  will  be  lost.  Some  of  these  prerequisites  are  not  necessary  when 
modeling  steady  flow,  but  are  critical  for  unsteady  turbulent  flow.  The  cases  in  this  chapter 
are  designed  to  ensure  that  these  prerequisites  are  met. 

Three  non-turbulent  cases  will  be  examined;  a  vortex  convecting  in  an  inviscid  atmo¬ 
sphere,  a  viscous  shock-tube,  and  a  2D  cylinder  in  crossflow.  As  described  in  Section  2.2, 
the  Roe  and  Steger- Warming  methods  are  included  in  Beggar  to  solve  the  right  hand  side 
term,  which  is  a  primary  factor  of  artificial  dissipation.  Each  of  the  three  cases  will  examine 
the  differences  between  the  Roe  and  Steger- Warming  discretization  schemes  to  approach  a 
recommendation  on  which  should  be  used  in  unsteady  turbulent  cases.  The  best  method 
will  be  used  in  the  three  turbulent  cases  in  Chapter  IV,  because  the  time  needed  to  run 
these  cases  for  both  schemes  would  exceed  the  time  available  to  finish  this  study. 

As  described  in  Section  2.2,  Beggar  uses  a  Newton  iteration  scheme  to  converge  each 
time  step.  The  shock-tube  and  2D  cylinder  cases  were  used  to  evaluate  various  numbers 
of  Newton  dt-iterations.  Determining  the  number  of  dt-iterations  is  an  important  aspect 
of  the  non-turbulent  cases.  Unlike  steady-state  calculations,  when  converging  each  time 
step  is  not  an  important  consideration,  each  time  step  must  be  converged  to  properly  sim¬ 
ulate  unsteady  flow.  For  unsteady  Beggar  computation,  the  user  specifies  both  the  number 
of  Newton  dt-iterations  and  the  maximum  number  of  Gauss-Seidel  inner  iterations.  Con¬ 
vergence  criteria  may  also  be  set  on  the  Gauss-Seidel  inner  iterations.  Unless  there  is  a 
good  deal  of  variation,  for  example  initially  imposing  the  wall  boundary  conditions,  the 
maximum  number  of  Gauss-Seidel  inner  iterations  is  rarely  used.  However,  the  specified 
number  of  Newton  dt-iterations  is  always  used.  Each  Newton  iteration  is  computationally 
expensive.  The  amount  of  accuracy  gained  from  each  additional  Newton  dt-iteration  drops 
logarithmically. 

All  three  of  the  non-turbulent  cases  are  two-dimensional.  However,  Beggar  is  a  three- 
dimensional  code.  Therefore,  all  three  grids  actually  have  dimensions  of  nxmx 3.  To  force 
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a  2D  computation,  tangent  boundary  conditions  were  applied  to  both  outside  layers,  which 
means  that  the  flow  equations  were  only  solved  at  the  middle  layer  of  nodes. 


3.1  Invisid  Convecting  Vortex 

3.1.1  Case  Presentation.  The  first  case  examined  in  the  process  of  validating  the 
turbulence  model  is  a  vortex  convecting  in  an  inviscid  atmosphere.  The  inviscid  atmosphere 
is  modeled  using  the  Euler  equations  which  do  not  include  any  viscous  terms.  Therefore, 
any  dissipation  can  be  attributed  to  numerical  fluxes,  and  is  artificial.  For  steady  state 
solutions,  this  artificial  dissipation  may  be  desirable.  The  increased  dissipation  may  damp 
out  undesirable  unsteady  effects.  However,  time-accurate  solutions  are  dependent  on  these 
unsteady  effects,  and  too  much  artificial  dissipation  could  corrupt  the  simulation.  The 
inviscid  convecting  vortex  is  a  quick  simulation  to  run  to  check  the  numerical  dissipation  in 
each  component  of  the  code. 

A  traversing  inviscid  vortex  case  was  simulated  several  times  while  using  different 
solver  option  combinations  each  time.  The  following  solver  option  combinations  were  used: 


1.  2nd  order  Roe  right  hand  side  with  1st  order  backwards  Euler  time  step 

2.  2nd  order  Steger- Warming  right  hand  side  with  1st  order  backwards  Euler  time  step 

3.  2nd  order  Roe  right  hand  side  with  2nd  order  three  point  backwards  time  step 

4.  2nd  order  Steger- Warming  right  hand  side  with  2nd  order  three  point  backwards  time 
step 

The  first  order  backwards  Euler  time  scheme  is  less  accurate  than  the  second  order  three 
point  backwards  and  is  expected  to  have  more  numerical  dissipation.  Similarly,  the  Roe 
method  is  more  accurate  than  Steger- Warming  and  should  provide  less  numerical  dissipa¬ 
tion. 

Assuming  a  unit  vortex  core  radius,  the  non-dimensional  equations  for  generating  the 
initial  velocities  and  pressure  are(24): 

u  =  uoo  -  T^exp  Q(!  “  r2)^J  aV  (52) 
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and  R  is  the  distance  from  the  center  of  the  vortex,  is  the  speed  of  sound  and  Ax  and 
Ay  are  the  x  and  y  positions  with  respect  to  the  center  of  the  vortex,  respectively.  The 
strength  of  the  vortex  is  given  by  T.  The  initial  conditions  used  are  given  in  Table  2. 


Table  2:  Inviscid  convecting  vortex  initial  conditions  -  non-dimensional  values 


Property 

Value 

Aioo 

0.5 

Re 

2.0xl0t> 

r 

5.0 

Coo 

1.0 

Poo 

1.0 

Too 

1.4 

Solver  Options: 

2nd  Order  Spatial  Discretization 

Steger- Warming  Jacobians 

The  effect  of  a  vortex  is  negligible  at  a  distance  of  5  times  the  core  radius  from  the 
center  of  the  vortex.  Therefore,  with  a  vortex  core  radius  of  1.0,  the  effects  of  the  vortex 
will  be  negligible  at  a  distance  of  5.0  from  the  center.  From  this  knowledge,  grids  with  equal 
length  sides  of  10  vortex  core  radii  were  created.  All  cases  were  run  on  both  a  coarse  grid 
and  a  fine  grid.  The  coarse  grid  has  dimensions  of  81x81  while  the  fine  grid  has  dimensions 
of  161x161.  The  general  spacing  of  the  coarse  grid  with  respect  to  the  vortex  can  be  seen  in 
Figure  3.  However,  since  Beggar  is  a  cell  centered  code,  all  of  its  output  is  also  cell  centered. 
Therefore,  the  gridlines  seen  in  Figure  3  intersect  at  cell  centers  rather  than  gridpoints.  This 
gives  the  illusion  that  the  initial  center  of  the  vortex  lies  in  a  cell  center,  when  actually  the 
initial  center  of  the  vortex  lies  at  the  gridpoint  with  coordinates  (5. 0,5.0).  In  hindsight,  the 
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Figure  3:  Initial  vortex  in  the  coarse  81x81  grid 


vortex  center  should  have  been  placed  at  a  cell  center.  That  way  Beggar  would  have  had 
the  capability  of  exactly  predicting  the  y-location  of  the  center  of  the  vortex.  With  the 
center  of  the  vortex  being  placed  at  a  node  and  Beggar  being  a  cell-centered  code,  the  best 
Beggar  can  do  is  to  predict  the  center  of  the  vortex  77  cell  from  its  true  location. 

As  mentioned  previously,  the  vortex  was  initially  centered  at  the  middle  of  the  grid. 
The  vortex  then  convected  in  the  x  direction  at  the  set  Mach  number.  Periodic  boundary 
conditions  were  applied  on  all  boundaries.  Therefore,  when  the  vortex  passes  through  one 
boundary  and  instantaneously  appears  on  the  opposite  boundary.  A  cycle  is  complete  when 
the  vortex  returns  to  its  original  position.  Two  time  steps  were  investigated.  The  larger  and 
smaller  time  steps  have  non-dimensional  values  of  0.04  and  0.02,  respectively.  With  these 
time  step  sizes,  500  and  1000  time  steps  respectively  are  necessary  for  the  vortex  to  traverse 
the  length  of  the  grid  and  return  to  its  original  position.  Each  time  step  was  converged 
using  three  Newton  dt-iterations.  Each  Newton  dt-iteration  was  converged  using  at  most 


28 


60  symmetric  Gauss-Seidel  inner  iterations.  If  the  inner  iterations  reached  the  convergence 
criteria,  then  not  all  60  iterations  were  calculated. 

Pressure  was  used  to  indirectly  measure  the  errors  in  the  solution.  As  can  be  seen  in 
Equations  54  and  55,  the  pressure  is  minimum  at  the  center  of  the  vortex.  After  each  time 
step,  the  minimum  pressure  was  located.  The  magnitude  and  location  were  both  stored. 
Therefore,  there  are  two  items  to  measure  error: 

1.  Change  in  minimum  pressure  from  initial  value 

2.  Deviation  of  the  minimum  pressure  from  theoretical  location  of  vortex  center 

3.1.2  Results.  The  results  for  the  inviscid  vortex  case  are  shown  in  Figures  4  to 
18.  Figures  4  to  8  show  data  at  an  exact  point  in  time.  Figures  9  to  18  show  continuous 
data  from  zero  to  six  cycles. 

In  Figures  4  to  8,  only  the  contours  and  the  grid  outline  appear;  the  grid  has  been 
removed  to  give  a  better  view  of  the  contours.  Figure  4  shows  the  initial  pressure  contours 
and  the  direction  of  vortex  convection.  Figures  5  to  8  show  pressure  contours  after  six  cycles 
for  each  of  the  solver  combinations.  An  exact  solution  would  result  in  an  unchanged  vortex 
after  each  cycle.  Artificial  dissipation  can  be  qualitatively  seen  by  the  expansion  of  the  inner 
ring  of  the  pressure  contour  plots.  Other  numerical  errors  can  be  seen  in  the  distortion  of 
the  shape  of  the  vortex. 

Figure  5  shows  the  1st  order  temporal  -  Steger- Warming  RHS  combination.  It  has  a 
good  deal  of  dissipation  as  the  innermost  pressure  contour  ring  has  expanded  a  relatively 
large  amount.  The  shape  of  the  vortex  has  also  been  moderately  distorted.  The  vortex 
shape  resembles  more  of  an  ellipse  than  a  circle. 

Figure  6  shows  the  2nd  order  temporal  -  Steger- Warming  RHS  combination.  It  has 
an  innermost  pressure  contour  almost  identical  to  the  initial  vortex.  The  circular  shape  is 
maintained  throughout  most  of  the  core  of  the  vortex.  However,  the  outermost  rings  of  the 
vortex  are  highly  distorted. 

Figure  7  shows  the  1st  order  temporal  -  Roe  RHS  combination.  Similar  to  Figure 
5,  it  has  noticable  artificial  dissipation  as  seen  in  the  expansion  of  the  innermost  pressure 
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contour.  Unlike  Figure  5,  however,  Roe  did  not  distort  the  general  shape  of  the  vortex.  It 
still  much  resembles  a  circle  until  the  outermost  contour. 

Figure  8  shows  the  2nd  order  temporal  -  Roe  RHS  combination.  This  combination 
has  compacted  the  whole  vortex,  reducing  the  core  radius  and  strength  of  the  vortex.  The 
circular  shape  of  most  of  the  inner  contours  was  maintained.  The  outermost  contours  are 
distorted,  no  longer  resembling  a  circle  in  the  least.  Furthermore,  the  position  of  the  vortex  is 
noticably  different.  The  vortex  appears  to  be  centered  along  the  abscissa,  but  has  traversed 
about  one  vortex  core  radius  along  the  ordinate. 
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Figure  4:  Vortex  -  Initial  Pressure  Contours 
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Figure  5:  Vortex  -  Steger-Warming  -  1st  Order  Temporal  -  Pressure  Contours  After  6 
Cycles 


Figure  6:  Vortex  -  Steger-Warming  -  2nd  Order  Temporal  -  Pressure  Contours  After  6 
Cycles 
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Vortex  -  Roe  -  1st  Order  Temporal  -  Pressure  Contours  After  6  Cycles 
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Vortex  -  Roe  -  2nd  Order  Temporal  -  Pressure  Contours  After  6  Cycles 
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Figures  9  to  18  show  continuous  data  from  the  initial  vortex  to  when  the  vortex  has 
traversed  six  cycles.  At  each  time  step,  the  minimum  pressure  and  its  x  and  y  coordinates 
were  found.  Figures  9  and  10  show  the  minimum  pressure  as  the  vortex  traverses.  Recalling 
that  the  minimum  pressure  should  be  located  at  the  center  of  the  vortex,  the  x  and  y 
coordinates  of  the  minimum  pressure  at  each  time  step  were  used  as  the  vortex  center  to 
find  the  Ax  and  Ay  of  the  vortex  from  its  theoretical  location.  Figures  11  to  18  show 
various  deviations  as  the  vortex  traverses.  Because  the  true  center  of  the  vortex  passes 
through  grid-points  lying  on  the  y  =  5.0  axis  and  Beggar  is  a  cell-centered  instead  of  node- 
centered  scheme,  the  best  Ay  obtainable  is  either  Ay  =  +0.5  cells  or  Ay  =  —0.5  cells.  The 
computed  x  values  as  well  lie  at  discrete  points  whereas  the  theoretical  vortex  center  will  lie 
between  grid-points.  Therefore,  the  best  Ax  obtainable  will  range  from  —0.5  <  Ax  <  0.5 
cells,  depending  on  the  theoretical  position  of  the  vortex. 

Figure  9  compares  Roe  and  Steger- Warming,  1st  and  2nd  order  combinations  on  the 
coarse  grid.  The  initial  slope  of  the  curves  are  greater  for  both  1st  order  cases,  indicating 
more  numerical  dissipation.  This  difference  in  slope  shows  the  necessity  of  a  second  order 
time  discretization.  Within  the  first  two  cycles,  the  error  of  the  1st  order  time  discretization 
on  the  coarse  mesh  was  about  two  times  greater  than  the  2nd  order  time  discretization  with 
the  Steger- Warming  method  and  about  four  times  greater  with  the  Roe  method.  Comparing 
the  different  RHS  methods,  Roe  has  a  smaller  initial  slope  than  Steger- Warming  for  both  the 
1st  and  2nd  order  temporal  discretization  schemes.  Similar  results  were  obtained  with  the 
fine  grid,  shown  in  Figure  10,  except  at  about  five  cycles  Roe  with  the  2nd  order  temporal 
scheme  has  a  drastic  jump  in  minimum  pressure  and  apparently  becomes  unstable.  Similar 
results  were  observed  with  the  halved  time  step.  All  plots  can  be  seen  in  Appendix  A.l. 

Figures  11  and  12  shows  that  the  vortex  center  moved  -0.75  cells  in  the  x  direction 
while  staying  as  on  track  as  numerically  possible  in  the  y  direction  for  the  1st  order  temporal 
-  Roe  RHS  combination. 

Figures  13  and  14  shows  the  2nd  order  temporal  -  Roe  RHS  combination.  The  vortex 
starts  to  traverse  a  little  slower  than  theory  until  the  fifth  cycle.  Then  the  vortex  suddenly 
speeds  up  and  the  Ax  goes  from  one  cell  behind  to  3  cells  ahead  of  the  theoretical  location. 
Similarly,  the  y  position  is  as  true  as  possible  until  the  sixth  cycle.  Then  the  vortex  starts 
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traversing  in  the  negative  y  direction.  The  calculated  vortex  center  moved  more  than  six 
cells  in  the  y  direction  between  the  fifth  and  sixth  cycle.  This  shows  that  the  2nd  order 
temporal  -  Roe  RHS  combination  has  inherent  instabilities.  Roe  treats  2 D  problems  as  two 
ID  problems,  resulting  in  maximum  error  45°  from  these  axes(41).  This  combined  with  the 
second  order  temporal  discretization  could  be  causing  the  instability. 

Figures  15  and  16  shows  that  the  vortex  center  moved  an  average  of  -0.5  cells  in  the 
x  direction  for  the  1st  order  temporal  -  Steger- Warming  RHS  combination.  The  calculated 
y  position  was  predicted  perfectly  until  the  sixth  cycle,  when  it  moved  one  cell  down  from 
the  theoretical  position. 

Figures  17  and  18  shows  that  the  vortex  center  deviated  about  one  cell  in  the  x 
direction  in  five  cycles  before  starting  to  return  to  its  true  position  for  the  2nd  order  temporal 
-  Steger- Warming  RHS  combination.  The  calculated  y  position  was  predicted  perfectly 
except  near  the  completion  of  the  fifth  cycle  when  it  moved  one  cell  down  from  the  theoretical 
position. 

The  fine  grid  and  the  halved  time-step  simulations  gave  similar  results,  which  can  be 
seen  in  Appendix  A.l. 
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Figure  11:  Cycle  vs  Ax  -  Roe  -  Coarse  Grid  -  1st  Order  Temporal 
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Figure  13:  Cycle  vs  Ax  -  Roe  -  Coarse  Grid  -  2nd  Order  Temporal 
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Figure  15:  Cycle  vs  Ax  -  Steger- Warming  -  Coarse  -  1st  Order  Temporal 
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Figure  17:  Cycle  vs  Ax  -  Steger- Warming  -  Coarse  Grid  -  2nd  Order  Temporal 
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Figure  18:  Cycle  vs  Ay  -  Steger- Warming  -  Coarse  Grid  -  2nd  Order  Temporal 


39 


3. 2  Shock-  Tube 


3.2.1  Case  Presentation.  The  second  case  in  the  non-turbulent  study  is  the  shock- 
tube  or  ’Riemann’  problem(39).  The  theoretical  data  for  the  shock-tube  is  an  exact  solution 
of  the  Euler  equations  that  can  be  obtained  from  any  compressible  flow  textbook.  This  case 
is  interesting  because  it  contains  a  moving  shock,  a  growing  expansion  fan,  and  a  contact 
discontinuity.  To  slightly  increase  the  difficulty  of  the  problem,  viscous  laminar  Navier- 
Stokes  equations  were  solved  in  the  simulation.  This  was  accomplished  by  having  a  two 
dimensional  grid  with  a  wall  boundary  condition  on  the  top  and  bottom.  Figure  19  shows 
the  viscous  shock-tube  grid,  which  has  dimensions  of  101x81.  The  viscous  spacing  off  the 
axial  wall  can  also  be  seen  in  Figure  19,  which  has  an  initial  wall  spacing  of  0.09%  of  the 
total  height  of  the  grid.  The  center  of  the  shock  tube  along  the  axial  direction  is  far  enough 
from  the  walls  that  viscous  effects  will  be  minimal  and  will  be  nearly  inviscid.  Data  is 
gathered  from  this  location  as  it  will  best  match  the  inviscid  theoretical  data. 


Figure  19:  2D  Viscous  Shock- Tube  Grid 

The  shock-tube  is  simulated  by  dividing  the  long  rectangular  grid  shown  in  Figure 
19  into  left  and  right  regions.  The  left  and  right  regions  are  then  initialized  using  different 
values  for  density,  pressure,  and/or  temperature.  Both  regions  are  initialized  to  have  zero 
initial  velocity.  This  imitates  a  shock-tube  directly  after  the  diaphragm  has  been  broken. 
The  non-dimensional  values  used  in  this  case  are  shown  in  Table  3. 
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Table  3:  Shock-tube  example  initial  conditions  -  non-dimensional  values 


Property 

Left  Region 

Right  Region 

P 

1.0 

0.1 

P 

1.0 

0.1 

Mach 

0.0 

0.0 

Solver  Options: 

2nd  Order  Spatial  Discretization 
Steger-Warming  Jacobians 

The  shock-tube  case  was  used  to  perform  the  three  following  items: 

1.  Investigate  differences  between  the  Roe  and  Steger- Warming  discritization  schemes 

2.  Determine  the  optimal  number  of  Newton  dt-iterations  to  use  in  future  calculations 

To  test  that  the  solution  is  independent  of  the  time  step,  two  time  steps  will  be  used. 
The  larger  and  smaller  non-dimensional  time  steps  are  0.00125  and  0.000625,  respectively. 
All  data  was  recorded  at  a  non-dimensional  time  of  0.2  which  corresponds  to  160  and  320 
time  steps,  respectively. 

The  most  important  aspect  to  this  case  is  the  Newton  dt-iteration  test.  This  case 
was  run  with  one,  three,  five,  ten,  twenty,  and  fifty  Newton  dt-iterations.  The  fifty  Newton 
dt-iterations  runs  were  used  as  reference  solutions  to  compute  errors. 

3.2.2  Results.  The  results  for  the  shock-tube  case  are  shown  in  Figures  20  to  24. 

Figure  20  shows  a  density  plot  at  non-dimensional  time  0.2  along  the  centerline  of  the 
shock  tube.  Both  the  1st  and  the  2nd  order  time  discretization  methods  accurately  predict 
the  location  of  the  expansion  fan,  contact  discontinuity  and  the  shock.  Roe  does  a  better  job 
than  Steger-Warming  at  capturing  the  shock  and  contact  discontinuities,  and  more  sharply 
predicts  the  beginning  and  end  of  the  expansion  fan.  Figure  21  shows  a  zoomed  in  view  of 
the  trailing  edge  of  the  expansion  fan. 

Figure  22  shows  that  halving  the  time  step  and  doubling  the  number  of  time  iterations 
will  result  in  the  same  solution,  proving  that  the  converged  time  step  size  has  been  used. 
Figure  23  compares  a  1  Newton  dt-iteration  per  time  step  simulation  to  a  50  Newton  dt- 
iterations  per  time  step  simulation.  For  this  simple  problem,  1  Newton  dt-iteration  results 
in  as  good  of  a  solution  as  50  Newton  dt-iterations. 
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Figures  24  shows  the  results  of  the  Newton  dt-iteration  study.  For  this  study,  the 
50  Newton  dt-iterations  per  time  step  simulation  was  used  as  a  reference  to  calculate  the 
percent  errors  of  the  other  simulations  using  an  L ^  norm  calculation.  Instead  of  being  an 
error  with  respect  to  theory  or  experimental  data,  this  error  is  Newton  convergence  error. 
Furthermore,  as  can  be  seen  in  Figure  23,  both  the  1  and  50  Newton  dt-iterations  per  time 
step  simulations  gave  almost  the  same  answer.  Therefore,  the  error  plots  presented  here 
are  used  only  to  determine  Newton  dt-iteration  convergence,  not  to  determine  if  a  certain 
number  of  Newton  dt-iterations  is  within  an  error  bounds.  All  of  the  curves  in  Figure  24 
show  the  same  trends;  they  all  plateau  at  5  Newton  dt-iterations.  Error  data  gathered 
from  pressure  was  plotted  and  gave  similar  results.  It  is  concluded  that  more  than  5  dt- 
iterations  is  superfluous.  The  error  curve  breaks  between  3  and  5  Newton  dt-iterations. 
Before  3  Newton  dt-iterations  the  the  curve  is  steep  enough  that  using  less  than  3  Newton 
dt-iterations  per  time  step  may  result  in  poorly  converged  time  steps.  The  conclusion  is 
reached  that  between  3  and  5  Newton  dt-iterations  should  be  calculated  per  time  iteration. 


Figure  20:  Shock-Tube  Density  at  Non-Dimensional  Time  t  =  0.2 
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Figure  21:  Zoomed  in  Image  of  Figure  20  at  Trailing  Edge  of  Expansion  Fan 


Figure  22:  Shock-Tube  Density  at  Non-Dimensional  Time  t  =  0.2  -  Halved  Time  Step 
comparison 
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Figure  23:  Shock-Tube  Density  at  Non-Dimensional  Time  t  =  0.2  -  Newton  comparison 


Figure  24:  Maximum  Error  in  Density  vs  Number  of  Newton  dt-iterations 
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3.3  2D  Vortex  Shedding  Cylinder 


3.3.1  Case  Presentation.  The  shock-tube  provided  insight  into  the  ability  of  the 
code  to  accurately  solve  time  resolved  problems,  but  it  is  a  very  simple  case.  In  order  to 
more  fully  validate  the  code,  more  complex  problems  need  to  be  considered.  Shedding  of 
vortices  from  a  laminar  cylinder  in  crossflow  is  a  more  complex  flow  over  larger  time  scales. 
Experimental  data  shows(43)  that  for  Rer>  of  150  the  Strouhal  number  (St)  has  a  value 
slightly  less  than: 

fD 

St  =  -jj-  ~  0.2  (56) 

The  vortex  shedding  period,  A,  may  then  be  put  in  terms  of  Strouhal  number  via: 


A  = 


U 

1 


D 

St 


(57) 


A  good  visual  check  to  see  if  the  proper  Strouhal  number  was  calculated  is  to  see  if  the 
distance  between  vorticies  is  about  equal  to  five  cylinder  diameters.  Furthermore,  for  Re  d  < 
180,  viscous  effects  quickly  damp  any  flow  velocity  in  the  axial  direction,  reducing  the 
problem  to  two  dimensions.  This  makes  for  a  relatively  quick  simulation  with  valuable 
results.  The  conditions  chosen  for  this  simulation  are  shown  in  Table  4.  A  Newton  dt- 
iteration  convergence  study  was  performed  for  this  case.  The  following  number  of  dt- 
iterations  were  examined:  1,  3,  5,  10,  and  20.  Each  Newton  dt-iteration  solved  60  symmetric 
Gauss-Seidel  inner  iterations  unless  an  inner  convergence  tolerance  of  1  •  10~8  was  met. 


Table  4:  Vortex  shedding  laminar  cylinder  example  initial  conditions 


Property 

Value 

M 

0.2 

Ren 

150 

Tref 

500°i? 

At 

9.12  •  10“5  s 

Solver  options: 

Steger- Warming  Jacobians 

2nd  Order  Spatial  Differencing 

The  time  step  used  resulted  in  about  250  time  steps  between  consecutively  shedding 
vortices.  This  number  provides  enough  data  to  statistically  measure  St.  The  grid  used 
has  dimensions  of  401x201  and  is  shown  in  Figure  25.  The  cylinder  has  a  non-dimensional 
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Figure  25:  2D  Vortex  Shedding  Laminar  Cylinder  grid. 


diameter  of  1.0.  The  grid  extends  out  100  diameters  with  an  initial  spacing  of  0.001  diame¬ 
ters.  Similar  to  the  laminar  shock-tube,  two  different  temporal  discretization  methods  were 
evaluated,  a  1st  and  2nd  order  scheme.  This  case  was  run  for  10,000  iterations.  The  last 
4096  were  statistically  analyzed. 
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3.3.2  Results.  The  results  for  the  2D  vortex  shedding  cylinder  case  are  shown  in 
Figures  26  to  31. 

Figure  26  shows  an  example  contour  rnach  plot  of  the  cylinder.  The  plot  was  made  with 
data  from  using  Roe  RHS  with  2nd  order  temporal  discretization.  All  scheme  combinations 
provided  similar  figures.  Shedding  vortices  are  clearly  visible  and  a  visual  measurement 
of  Strouhal  number  (Period  ~  5  diameters)  looks  reasonable.  The  data  was  statistically 
analyzed  to  obtain  more  precise  measurements. 

Lift,  drag,  and  time  were  recorded  at  each  time  step.  An  example  plot  of  lift  vs  time  is 
shown  in  Figure  27.  Fast  Fourier  Transforms  (FFT)  were  performed  on  the  lift  data  to  find 
the  frequency  of  oscillation  via  a  power  spectrum  density  (PSD).  The  Strouhal  number  is 
determined  from  the  largest  peak  on  the  PSD  plot.  The  frequency  output  from  the  FFT  was 
non-dimensionalized  by  the  cylinder  diameter  and  freestream  velocity  to  form  a  Strouhal 
number.  Figures  28  to  30  are  a  few  of  these  PSD  plots. 

Figure  28  compares  the  temporal  order  discretization  methods.  Increased  dissipation 
from  numerical  error  stretches  the  shedding  period.  Stretching  the  shedding  period  results 
in  a  shortened  frequency,  which  in  turn  predicts  a  small  Strouhal  number.  The  first  order 
temporal  discretization  method  slightly  under-predicts  Strouhal  number  compared  to  the 
second  order  discretization  and,  therefore,  is  showing  more  numerical  dissipation  than  the 
second  order  temporal  discretization  method. 

Figure  29  compares  the  different  RHS  methods.  Using  the  same  rationale  as  with  the 
temporal  comparison  in  the  previous  paragraph,  it  is  apparent  that  Steger- Warming  has 
more  numerical  dissipation  than  Roe. 

Figure  30  shows  a  Newton  dt-iteration  comparison  for  the  Roe  RHS  -  1st  order  tem¬ 
poral  solver  option  combination.  Three  dt-iterations  per  time  step  and  up  all  predicted 
the  same  Strouhal  number  slightly  less  than  0.2.  However,  the  1  dt-iteration  per  time  step 
simulation  case  predicted  a  Strouhal  number  lower  than  all  the  other  cases.  Similar  results 
were  found  with  other  solver  option  combinations  as  are  shown  in  Appendix  A. 2. 

The  maximum  lift  was  recorded  for  each  shedding  vortex,  and  the  average  value  was 
found.  Figure  31  shows  average  peak  drag  with  respect  to  Newton  dt-iterations.  In  general, 
one  dt-iteration  is  far  from  a  converged  solution  and  three  dt-iterations  is  well  converged. 
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The  solution  is  always  converged  for  five  dt-iterations  or  greater.  Similar  plots  were  made 
with  other  properties  as  shown  in  Appendix  A. 2.  All  results  showed  the  same  trends. 
Therefore,  the  author  recommends  using  between  three  and  five  dt-iterations  to  properly 
converge  each  time  iterations  for  unsteady  turbulent  flow. 


Figure  26:  Mach  contours  of  vortex  shedding  laminar  cylinder 
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Figure  27: 


Figure  28:  PSD  of  St  from 


Time  (sec) 


Lift  vs  Time 


-  Temporal  Order  Comparison 


PSD  Cx 


Figure  29:  PSD  of  St  from  Lift  -  RHS  Comparison 


Figure  30:  PSD  of  St  from  Lift  -  Roe  -  1st  Order  Temporal 


Figure  31: 


Average  Peak  Lift  vs  dt-iteration 


51 


3.4  Conclusions 

To  properly  simulate  unsteady  turbulent  flow,  it  is  important  to  reduce  artificial  dis¬ 
sipation  and  accurately  calculate  each  time  step.  The  inviscid  convecting  vortex  provided 
useful  data  about  the  contribution  of  each  solver  option  to  artificial  dissipation,  while  the 
shock-tube  and  2D  cylinder  were  most  beneficial  in  the  area  of  time  step  convergence.  Re¬ 
viewing  the  non-turbulent  cases,  the  conclusions  are  drawn  that  the  following  solver  options 
should  be  implemented  in  the  Beggar  code  to  simulate  unsteady  turbulent  flow: 

1.  A  second  order  temporal  discretization 

2.  Between  three  and  five  Newton  dt-iterations 

The  Roe  method  ultimately  was  used  to  compute  the  RHS  for  the  cases  in  Chapter  4  because 
it  has  been  shown  to  have  less  numerical  dissipation.  However,  more  investigation  should 
be  done  on  the  instabilities  it  showed  in  the  inviscid  vortex  case. 
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IV.  TURBULENT  STUDY 


With  the  Beggar  code  validated  for  unsteady  flows,  the  turbulence  models  may  now  be 
inspected.  Many  of  the  turbulence  models  used  today  were  developed  for  steady-state 
RANS  flow  solvers.  They  include  significant  numerical  dissipation  to  wash  out  oscillations 
and  rapidly  converge  to  a  steady  state  solution.  This  numerical  dissipation  will  allow  a 
steady-state  problem  to  converge  much  quicker,  but  it  may  destroy  the  ability  of  the  code 
to  predict  unsteady  flow.  Therefore,  unsteady  cases  must  be  simulated  in  order  to  validate 
a  turbulence  model  for  unsteady  flow.  The  unsteady  turbulent  validation  cases  selected  for 
this  effort  are  a  2D  oscillating  airfoil,  a  3D  turbulent  shedding  cylinder,  and  a  turbulent 
cavity. 

4-1  2D  Oscillating  Airfoil 

4-1.1  Case  Presentation.  A  simple  turbulent  case  is  that  of  an  oscillating  (pitch¬ 
ing)  2D  airfoil.  The  time  dependance  of  the  angle  of  attack,  a,  is  given  by 

a  =  a0  +  4.2°sm(207rt)  (58) 

Two  af  s  were  investigated,  aQ  =  4.0°  and  a0  =  11.0°.  For  the  a0  =  4.0°  case,  the  a’s  are 
kept  small  so  that  the  flow  always  remains  attached  to  the  airfoil(13).  With  the  flow  always 
attached,  Baldwin-Lomax  would  be  expected  to  perform  as  well  as  Spalart-Allmaras  and 
DES.  For  the  a0  =  11.0°  case,  the  flow  separates  on  the  upper  surface  of  the  airfoil.  Here, 
Spalart-Allmaras  would  be  expected  to  outperform  Baldwin-Lomax. 

A  NACA  0015  airfoil  with  a  non-dimensional  chord  of  1.0  was  discretized  using  two 
overlapping  grids.  The  region  closest  to  the  airfoil  is  covered  by  a  C-grid  with  dimensions 
381x55.  This  airfoil  grid  extends  1  chord  away  from  the  airfoil  and  6  chords  behind  the 
airfoil.  The  initial  non-dimensional  spacing  off  the  airfoil  is  Ax  =  0.00001.  Figure  32  shows 
the  airfoil  grid  and  Figure  33  shows  the  airfoil  zoomed  in  at  the  trailing  edge.  To  oscillate 
the  airfoil,  the  airfoil  grid  is  rotated  inside  a  larger  stationary  cartesian  grid.  This  grid  has 
dimensions  141x81  with  lengths  of  31x20  chord  lengths.  Figure  34  shows  the  cartesian  grid 
with  cells  cut  out  where  the  airfoil  resides.  To  achieve  time  accurate  results,  five  Newton 
dt-iterations  were  run  every  time  iteration.  Each  Newton  dt-iteration  solved  80  symmetric 
Gauss-Seidel  inner  iterations  unless  an  inner  convergence  tolerance  of  1  •  10~8  is  met. 
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Figure  32:  NACA  0015  Grid 


Turbulence  models  to  be  examined  include  Baldwin-Lomax,  Spalart-Allmaras,  and  a 
DES  model.  The  input  conditions  are  shown  in  Table  5. 

Table  5:  2D  oscillating  airfoil  example  initial  conditions 


Property 

Value 

M 

0.29 

Rec 

1.95  •  10B 

At  (aQ  =  4.0°) 

2.24  •  10“4  s 

At  (a0  =  11.0°) 

3.125  •  10~5  s 

Solver  options: 

Steger- Warming  Jacobians 

2nd  Order  Spatial  Differencing 
Roe  Right  Hand  Side 

5  dt-iterations 

To  ensure  time  convergence  is  not  an  issue,  two  time  step  sizes  were  investigated  with 
the  a0  =  4.0°  case.  The  time  step  of  At  =  2.24  •  10-4  was  halved,  and  twice  as  many 
iterations  should  return  the  same  solution. 

The  experimental  data  was  gathered  from  reference  (27).  Pressure  transducers  were 
placed  at  four  stations  along  the  span  of  a  60  inch  semispan  NACA  0015  wing  with  a  12.0 
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Figure  33:  NACA  0015  Grid  -  Zoomed  in  on  Trailing  Edge 


Figure  34:  Cartesian  Grid 
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inch  chord  with  zero  twist.  The  main  station  located  at  mid-span  included  20  pressure 
transducers  along  the  chord.  The  three  off-center  stations  were  used  to  check  the  two- 
dimensionality  for  each  data  point.  Measurements  were  taken  240  times  per  oscillation. 
The  pressure  transducer  output  was  then  sent  through  a  500  Hz  low  pass  filter.  Because 
the  experimental  data  was  only  taken  at  20  points  along  the  chord  and  then  sent  through 
a  filter,  the  net  effect  of  the  experimental  setup  was  that  the  unsteady  data  are  highly 
filtered(23). 
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4-1.2  Results.  Results  for  the  2D  oscillating  airfoil  are  shown  in  Figures  35  through 
46.  Figures  35  through  40  apply  to  the  aQ  =  4.0°  simulations  while  Figures  41  through  46 
apply  to  the  a0  =  11.0°  simulations.  The  a0  =  4.0°  case  is  an  attached  flow  case(13) 
whereas  the  a0  =  11.0°  case  has  separation  near  the  trailing  edge  as  the  airfoil  pitches  up 
(upstroke)  and  no  separation  as  the  airfoil  pitches  down  (downstroke). 

Figures  35  through  37  show  lift,  drag,  and  pitching  moment  coefficients  as  a  function 
of  angle  of  attack.  A  circuit  is  formed  on  each  due  to  the  induced  angle  of  attack  from  the 
oscillation.  All  models  accurately  calculated  the  lift  coefficient  as  seen  in  Figure  35.  Figure 
36  shows  all  three  of  the  models  under-predicting  C'o  at  a  =  —0.2°.  This  phenomena  has 
been  seen  in  other  codes  for  the  Baldwin-Lomax  and  Spalart-Allmaras  models(13,  35,  22). 
All  of  the  models  match  each  other  well  for  Cl  and  Co-  Baldwin-Lomax  diverges  from 
Spalart-Allmaras  and  DES  in  the  downstroke  portion  of  the  Cm  curve  as  shown  in  Figure 
37.  To  investigate  this,  plots  of  u- velocity  were  made  near  the  trailing  edge  of  the  airfoil  at 
a  =  5.22°  on  the  downstroke.  Figure  38  shows  the  Baldwin-Lomax  prediction  and  Figure 
39  shows  the  Spalart-Allmaras  prediction.  Baldwin-Lomax  has  a  large  region  of  negative 
velocity,  indicating  that  it  is  predicting  separation  on  the  trailing  edge.  The  flow  remains 
attached  in  the  Spalart-Allmaras  prediction,  as  there  is  no  negative  velocity.  As  noted  by 
Wilcox(45),  algebraic  models  usually  predict  too  much  separation.  It  appears  this  may  be 
the  case. 

Simulations  with  a  halved  time  step  were  also  run  to  see  if  the  solution  was  time 
converged.  Figure  40  shows  the  Cm  predicted  with  the  Baldwin-Lomax  model  for  the 
halved  time  step  compared  with  the  whole  time  step.  There  is  some  variation  at  high  angles 
of  attack  on  the  downstroke  where  Baldwin-Lomax  predicted  separation.  Otherwise,  the 
plots  lie  on  top  of  each  other.  Furthermore,  Spalart-Allmaras  and  DES  predictions  for  lift, 
drag,  or  moment  coefficient  did  not  change  when  the  time  step  was  halved.  It  is  therefore 
assumed  that  the  time  step  has  been  converged. 

Figures  41  through  43  show  lift,  drag,  and  pitching  moment  coefficients  as  a  function 
of  angle  of  attack  for  the  aQ  =  11.0°  simulations.  At  this  angle  of  attack  and  Reynolds 
number,  separation  is  expected  on  the  upstroke  and  disappears  on  the  downstroke.  This  is 


57 


known  as  a  light  dynamic  stall  case.  This  is  a  good  case  to  run  to  test  turbulence  models, 
as  light  dynamic  stall  is  sensitive  to  the  onset  of  separation. 

A  common  trend  among  the  three  models  is  that  they  all  performed  poorly  for  the 
aQ  =  11.0°  case.  None  of  the  models  resemble  the  lift  coefficient  curve  in  Figure  41.  All  the 
models  predicted  a  rounded  end  at  a  =  6.8°  and  a  point  at  a  =  15.2°,  where  data  shows 
the  opposite.  There  are  severe  oscillations  in  the  Baldwin-Lomax  curve,  indicating  that  the 
airfoil  is  shedding  vorticies. 

(35),  (13:Ko),  and  (23)  obtained  different  results  for  this  case  with  the  Spalart- 
Allmaras  model.  They  found  that  the  Spalart-Allmaras  model  matched  well  with  experi¬ 
mental  data  for  lift  and  drag,  and  had  the  same  shape  of  curve  for  Cm.  Ko  is  the  only  one 
who  also  compared  Baldwin-Lomax.  Ko’s  version  of  Baldwin-Lomax  predicted  the  same 
plot  on  the  upstroke  of  Figure  41,  but  did  not  predict  shedding  on  the  downstoke. 

Figures  44,  45,  and  46  show  vorticity  predictions  at  the  trailing  edge  of  the  airfoil  for 
Baldwin-Lomax,  Spalart-Allmaras,  and  DES  models,  respectively.  As  seen  in  these  figures, 
large  vortical  structures  are  shedding  from  the  airfoil  in  the  Baldwin-Lomax  simulation, 
whereas  no  vorticies  are  being  shed  from  either  the  Spalart-Allmaras  or  DES  simulations. 
Examining  the  velocity  vectors  showed  that  Baldwin-Lomax  predicted  separation  a  tenth 
of  a  chord  length  sooner  than  Spalart-Allmaras  and  DES.  This  distance,  combined  with 
the  fact  that  the  airfoil  has  more  curvature  closer  to  the  leading  edge  and  therefore  a  more 
adverse  pressure  gradient,  could  be  the  reason  Baldwin-Lomax  is  predicting  vortex  shedding 
while  Spalart-Allmaras  and  DES  are  not.  Recalling  that  Baldwin-Lomax  also  predicted 
separation  for  the  aQ  =  4.0°  case,  a  good  assumption  is  that  Baldwin-Lomax  is  predicting 
separation  too  early.  Reasons  why  Beggar’s  Spalart-Allmaras  data  does  not  match  well 
with  experimental  data  or  data  gathered  from  other  versions  of  the  Spalart-Allmaras  model 
cannot  be  determined,  but  the  question  is  raised  as  to  whether  the  Spalart-Allmaras  model 
was  coded  into  Beggar  properly. 

The  clock  times  were  recorded  for  the  aQ  =  4.0°  cases.  The  total  time  was  divided  by 
the  total  number  of  iterations  to  come  up  with  the  times  per  iteration  as  shown  in  Table  6. 
Because  Spalart-Allmaras  and  DES  have  to  solve  an  extra  differential  equation,  all  things 
being  equal  it  is  expected  that  those  simulations  would  take  longer  than  Baldwin-Lomax. 
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Table  6: 


Average  computation  time  for  various  turbulence  models  for  the  oscillating  airfoil 


case 


Model 

Time  per  Iteration 
Normalized  by  Baldwin-Lomax 

Baldwin-Lomax 

1.0 

Spalart-Allmaras 

1.27 

DES 

1.16 

Angle  of  Attack 

Figure  35:  a0  =  4.0°  -  Lift  Coefficient  vs  a 

However,  since  Baldwin-Lomax  predicted  shedding  vorticies  and  the  Spalart-Allmaras  model 
did  not,  the  Baldwin-Lomax  model  had  more  gradients  to  solve  for  and  the  number  of  inner 
iterations  needed  to  converge  each  Newton  dt-iteration  was  much  higher. 
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Figure  36:  aQ  =  4.0°  -  Drag  Coefficient  vs  a 
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Figure  37:  aD  =  4.0°  -  Pitching  Moment  Coefficient  vs  a 
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Figure  38: 


u  Velocity  -  a  =  5.22°  -  Downstroke-  Baldwin-Lomax 
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Figure  39:  u  Velocity  -  a  =  5.22°  -  Downstroke-  Spalart-Allmaras 
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Figure  40: 


Angle  of  Attack 

a0  =  4.0°  -  Pitching  Moment  Coefficient  vs  a  -  Half  Time  Comparison 
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Figure  41:  aQ  =  11.0°  -  Lift  Coefficient  vs  a 


Figure  42:  a0  =  11.0°  -  Drag  Coefficient  vs  a 
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Figure  43:  aQ  =  11.0°  -  Pitching  Moment  Coefficient  vs  a 


Figure  44:  Instantaneous  Vorticity  -  a0  =  11.0°  -  a  =  13.97°  Downstroke  -  Baldwin- 
Lomax 
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Figure  45:  Instantaneous  Vorticity  -  aQ  =  11.0°  -  a  =  13.97°  Downstroke  -  Spalart- 
Allmaras 
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Figure  46:  Instantaneous  Vorticity  -  aQ  =  11.0°  -  a  =  13.97°  Downstroke  -  DES 
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4-2  3D  Vortex  Shedding  Cylinder 


4-2.1  Case  Presentation.  A  vortex  shedding  turbulent  cylinder  is  significantly 
more  difficult  than  the  2D  case  examined  in  Chapter  III.  With  a  Re  >  180,  the  flow  can  no 
longer  be  approximated  as  2D.  A  vortex  may  be  shedding  from  the  top  of  the  cylinder  at 
one  point  along  the  span  on  the  cylinder  while  another  vortex  is  shedding  from  the  bottom 
of  the  cylinder  at  a  different  point.  Due  to  computational  and  time  constraints  and  some 
technical  difficulties,  only  a  coarse  mesh  study  was  performed.  With  the  coarse  mesh,  DES 
will  have  difficulty  performing  well.  With  the  larger  cells  in  the  wake,  only  the  largest  of 
the  eddies  will  be  captured  by  the  grid  and  DES  will  rely  on  the  sub-grid  model  to  account 
for  the  rest  of  the  eddies. 

The  dimensions  of  the  grid  are  101x73x51.  The  initial  spacing  off  the  wall  is  l.OxlO-'5 
which  gives  a  y+  ~  1.0.  The  cylinder  has  a  non-dimensional  diameter  of  1.0,  a  length  of 
10  chords,  and  the  grid  extends  50  chords  from  the  cylinder.  Initial  conditions  are  shown 
in  Table  7.  To  achieve  time  accurate  results,  five  Newton  dt-iterations  were  run  every  time 
step.  Each  Newton  dt-iteration  solved  80  symmetric  Gauss-Seidel  inner  iterations  unless  an 
inner  convergence  tolerance  of  1  •  10~8  was  met. 

Table  7:  Vortex  shedding  turbulent  cylinder  example  initial  conditions 


Property 

Value 

M 

0.2 

Ren 

8-10° 

At 

9.0  •  10"5  s 

Solver  options: 

Steger- Warming  Jacobians 

2nd  Order  Spatial  Differencing 
Roe  Right  Hand  Side 

5  dt-iterations 

4-2.2  Results.  Results  for  the  3D  vortex  shedding  cylinder  are  shown  in  Figures 
48  through  57  .  Pressure  coefficient  and  time  was  recorded  for  every  time  iteration  at  a 
point  on  top  of  the  cylinder.  This  data  was  passed  through  a  Fast  Fourier  Transform  to 
find  the  PSD’s  of  the  shedding  frequency.  The  shedding  frequency  was  non-dimensionalized 
which  results  in  Strouhal  number. 
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Figure  47  is  a  sample  pressure  coefficient  vs  time  plot  from  data  gathered  at  the  top 
of  the  cylinder.  Figure  48  shows  the  PSD  with  respect  to  St  computed  from  data  gathered 
at  the  top  of  the  cylinder.  Experiment(lO)  shows  for  Re £>  =  8  •  106, 

St  ss  0.306  -  0.308 

Figure  48  shows  both  the  Spalart-Allmaras  and  DES  models  underpredicted  the  Strouhal 
number  by  about  |.  Baldwin-Lomax  predicts  a  Strouhal  number  about  ^  the  value  of  the 
experimental  data.  The  smaller  turbulent  frequencies  in  the  Baldwin-Lomax  simulation  are 
overshadowing  the  vortex  shedding  frequency. 

Figures  49  to  51  show  instantaneous  plots  of  total  pressure  from  a  cross-sectional 
view  for  the  Baldwin-Lomax,  Spalart-Allmaras,  and  DES  models,  respectively.  In  each 
diagram  shedding  vorticies  are  present.  Spalart-Allmaras  and  DES  predict  periodic  shedding 
vorticies,  whereas  the  Baldwin-Lomax  model  predicts  randomly  shedding  vorticies. 

Figures  52  to  54  show  instantaneous  plots  of  vorticity  on  the  wake  surface  behind 
the  cylinder  and  the  wall  perpendicular  to  the  cylinder.  Three  dimensional  structures  are 
evident  in  the  Baldwin-Lomax  simulation,  as  seen  in  Figure  52.  However,  it  appears  that 
Spalart-Allmaras  and  DES  are  predicting  2D  vortex  shedding.  To  investigate  this,  cross 
sectional  plots  of  velocity  parallel  to  the  cylinder  were  made,  as  shown  Figures  55  to  57. 
Significant  velocity  parallel  to  the  cylinder  is  present  in  the  Baldwin-Lomax  simulation, 
Figure  55.  The  Spalart-Allmaras  and  DES  simulations  predict  almost  no  velocity  parallel 
to  the  cylinder,  as  seen  in  Figures  56  and  57  respectively. 

Table  8  shows  drag  coefficient  and  Strouhal  number  data  compared  to  experimental 
data  and  other  CFD  results.  The  average  Cd  matches  well  for  Baldwin-Lomax,  but  the 


Table  8:  Force  coefficient  and  Strouhal  number  data 


Model 

Average  Cc] 

St 

Baldwin-Lomax 

0.557 

0.11 

Spalart-Allmaras 

0.193 

.195 

DES 

0.189 

0.195 

Spalart-Allmaras  (25) 

0.226 

0.232 

DES  (25) 

0.298 

0.232 

Data  (10) 

0.505-0.540 

0.306-0.308 
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Time  (sec) 


Figure  47:  Cp  vs  Time 


Strouhal  number  was  underpredicted.  The  Spalart-Allmaras  and  DES  models  do  not  match 
the  average  Cd  or  St.  Using  the  same  grids  and  turbulence  model,  the  Spalart-Allmaras 
simulation  with  Beggar  predicted  a  Cd  15%  lower  than  predicted  with  the  NXAIR(25) 
code.  These  results  strengthen  the  belief  that  the  Spalart-Allmaras  model  was  not  properly 
written  in  the  Beggar  code. 

All  cases  were  run  on  a  1.8  GHz  Opteron  workstation  with  eight  gigabytes  of  memory. 
The  wall  clock  times  for  each  3D  cylinder  run  were  recorded.  The  average  time  per  iteration 
for  each  turbulence  model  is  shown  in  Table  9,  which  shows  DES  and  Spalart-Allmaras 
taking  about  the  same  time.  Baldwin-Lomax  took  about  1.5  times  longer  to  solve  each  time 
step.  As  described  previously,  Spalart-Allmaras  overdamped  the  unsteady  effects  which 
would  make  converging  each  time  step  simpler.  If  the  same  amount  of  unsteady  effects 
were  present  in  all  cases,  one  would  expect  Spalart-Allmaras  to  take  longer  to  solve  than 
Baldwin-Lomax. 

Table  9:  Average  computation  times  for  various  turbulence  models  for  the  3D  cylinder 
case  _ 


Model 

Time  per  Iteration 
Normalized  by  Baldwin-Lomax 

Baldwin-Lomax 

1.0 

Spalart-Allmaras 

0.62 

DES 

0.62 
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Strouhal  Number 

Figure  48:  PSD  of  St 


Figure  49:  Instantaneous  Total  Pressure  Plot  -  Baldwin-Lomax 
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Figure  50:  Instantaneous  Total  Pressure  Plot  -  Spalart-Allmaras 


Figure  51:  Instantaneous  Total  Pressure  Plot  -  DES 
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Figure  52:  Instantaneous  Vorticity  Plot  -  Baldwin-Lomax 


Figure  53:  Instantaneous  Vorticity  Plot  -  Spalart-Allmaras 
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Figure  54:  Instantaneous  Vorticity  Plot  -  DES 
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Figure  55:  Velocity  Parallel  to  the  Cylinder  -  Baldwin-Lomax 


73 


v_MA( 

0.200 


0.029 


-0.029 


-0.086 


Figure  56:  Velocity  Parallel  to  the  Cylinder  -  Spalart-Allmaras 
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Figure  57:  Velocity  Parallel  to  the  Cylinder  -  DES 


74 


4-3  Turbulent  Cavity 

4-3.1  Case  Presentation.  An  unsteady  simulation  of  the  Weapons  Internal  Car¬ 
riage  and  Separation  (WICS)  weapons  bay  was  also  investigated.  The  bay  has  the  physical 
dimensions  of  18x4x4  inches.  The  bay  was  placed  inside  grids  that  extend  15  inches  up¬ 
stream  of  the  bay  and  25  inches  downstream  of  the  bay.  This  matches  the  experimental 
setup(5).  A  grid  converge  study  was  performed  for  this  case.  Information  about  each  grid 
is  shown  in  Table  10  and  the  grids  are  shown  in  Figure  58. 


Table  10:  Turbulent  cavity  grid  specifications 


Grid 

Total  Points 

Bay  Grid 
Dimensions 

Bay  Grid 

Bay  Grid 

^  2/max 

Bay  Grid 
Ar 

^4max 

Fine 

1.9x10° 

133x137x61 

0.3  in 

0.1  in 

0.1  in 

Medium 

1.1x10° 

76x47x41 

0.6  in 

0.2  in 

0.2  in 

Coarse 

7.9xl05 

61x31x31 

0.75  in 

0.3  in 

0.3  in 

Initial  conditions  are  shown  in  Table  11.  To  achieve  time  accurate  results,  three  New¬ 
ton  dt-iterations  were  run  every  time  step.  Each  Newton  dt-iteration  solved  60  symmetric 
Gauss-Seidel  inner  iterations  unless  an  inner  convergence  tolerance  of  1  •  10~8  was  met. 

Table  11:  Turbulent  cavity  example  initial  conditions 


Property 

Value 

M 

0.95 

Re 

2.1  •  10s* 

At 

8.0  •  10~°  s 

Solver  options: 

Steger- Warming  Jacobians 

2nd  Order  Spatial  Differencing 
Roe  Right  Hand  Side 

3  dt-iterations 

Wind  tunnel  data  has  been  gathered  (5)  from  the  AEDC  four-foot  transonic  wind 
tunnel  using  two  transducers,  labeled  K16  and  K18.  Transducer  K-16  was  located  on  the 
centerline  of  the  bay  ceiling  0.275  inches  from  the  back  wall  while  K-18  was  located  on  the 
centerline  of  the  back  wall  0.725  inches  from  the  bay  opening.  CFD  data  was  gathered  at 
the  same  locations  and  compared  with  the  experimental  data. 
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Medium 


Fine 

Figure  58:  Bay  Grids 
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4-3.2  Results.  Results  for  the  turbulent  cavity  case  are  shown  in  Figures  59 
through  65.  Figures  59  through  62  compare  data  gathered  from  using  the  different  models 
to  the  experimental  data.  In  these  plots,  the  Spalart-Allmaras  model  was  far  from  matching 
the  experimental  data.  The  decible  reading  from  the  Spalart-Allmaras  model  in  Figure  59 
was  about  40  dB  less  than  the  experimental  data.  Only  a  few  frequencies  have  strong 
decibles  readings  with  the  Spalart-Allmaras  model.  This  indicates  that  only  the  largest 
fluctuations  are  strong  enough  to  propogate  all  the  way  to  the  bay  ceiling  and  all  the  other 
fluctuations  are  damped  out  by  the  model.  Baldwin-Lomax  and  DES  were  both  in  good 
agreement  with  the  experimental  data.  (37:Suhs)  also  found  Baldwin-Lomax  to  be  in  good 
agreement  with  the  experimental  data.  DES  does  a  better  job  at  matching  the  first  peak. 
This  simulation  was  also  run  by  (23:Nichols),  whose  data  shows  good  agreement  with  the 
Spalart-Allmaras  model  until  a  frequency  of  1000  Hz.  After  1000  Hz,  the  Spalart-Allmaras 
model  underpredicted  the  sound  level  by  about  20  decibles.  As  described  in  Section  2.7,  the 
only  difference  between  the  Spalart-Allmaras  and  DES  models  is  that  when  the  distance 
to  the  nearest  wall  gets  larger  than  0.65  times  the  largest  cell  dimension,  then  that  value 
is  used  in  place  of  the  distance  to  the  wall  to  solve  the  additional  PDE.  With  such  a  large 
difference  in  results  between  the  Spalart-Allmaras  and  DES  models,  it  is  likely  that  the 
distance  to  the  wall  is  not  being  calculated  properly. 

Figures  60  through  62  show  the  average  pressure  contour  along  the  centerline  of  the 
bay  ceiling  for  the  coarse,  medium,  and  fine  grids,  respectively.  DES  agrees  well  with 
the  experimental  data  when  solved  on  the  medium  grid,  but  predicts  almost  no  change 
in  pressure  coefficient  along  the  bay  ceiling  with  the  coarse  grid.  This  shows  that  the 
DES  model  is  more  sensitive  to  mesh  density  than  a  regular  RANS  model.  The  DES 
model  could  not  be  run  with  the  fine  grid  because  of  technical  difficulties.  Baldwin-Lomax 
performs  increasingly  well  as  the  mesh  density  is  increased.  For  the  fine  grid,  Baldwin- 
Lomax  data  matches  the  experimental  data  well.  For  each  of  the  coarse,  medium,  and  fine 
grid  simulations,  Spalart-Allmaras  predicts  zero  pressure  coefficient  along  the  centerline  of 
the  bay  ceiling  until  the  last  two  inches  of  the  bay.  Nichols(23)  shows  the  Spalart-Allmaras 
model  matching  the  pressure  coefficient  data  nicely,  making  the  case  stronger  that  the 
Spalart-Allmaras  model  in  Beggar  may  have  been  written  incorrectly  and  is  overdamping 
unsteady  effects. 
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Figure  59:  Sound  Pressure  Level  Spectrum  at  the  K16  Location  -  Medium  Grid 

Figures  63  through  65  show  instantaneous  vorticity  plots  on  the  centerplane  of  the 
bay  grid.  These  plots  give  a  better  visual  of  the  unsteady  effects,  or  lack  thereof  in  the 
Spalart-Allmaras  case.  Data  from  the  K18  point  and  all  other  grids  provided  similar  results 
and  can  be  seen  in  Appendix  B.l. 

All  runs  for  the  medium  case  were  performed  on  3  nodes  of  a  2.2  GHz  Opteron  cluster 
with  2  processors  per  node  and  four  gigabytes  of  memory  per  node.  The  wall  clock  times  for 
each  run  were  recorded.  The  average  time  per  iteration  for  each  turbulence  model  is  shown 
in  Table  12.  The  differences  in  time  are  not  significant  for  this  simulation.  The  additional 
time  needed  to  solve  the  extra  PDE  in  the  Spalart-Allmaras  model  balanced  the  additional 
time  needed  to  converge  the  Navier-Stokes  equations  with  the  Baldwin-Lomax  model  due 
to  the  extra  unsteady  effects  present  in  the  flow. 

Table  12:  Average  computation  times  for  various  turbulence  models  for  the  tubulent 
cavity  case 


Model 

Time  per  Iteration 
Normalized  by  Baldwin-Lomax 

Baldwin-Lomax 

1.0 

Spalart-Allmaras 

1.08 

DES 

0.99 
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Figure  60:  Average  Pressure  Coefficient  on  the  Cavity  Ceiling  -  Coarse  Grid 


Figure  61:  Average  Pressure  Coefficient  on  the  Cavity  Ceiling  -  Medium  Grid 
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Figure  62:  Average  Pressure  Coefficient  on  the  Cavity  Ceiling  -  Fine  Grid 
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Figure  63:  Vorticity  -  Baldwin-Lomax  -  Medium  Cavity  Grid 


80 


Figure  64:  Vorticity  -  Spalart-Allmaras  -  Medium  Cavity  Grid 
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Figure  65:  Vorticity  -  DES  -  Medium  Cavity  Grid 
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4-4  Conclusions 

Considering  all  three  cases  in  the  turbulent  study,  it  is  determined  that  the  version 
of  the  Spalart-Allnraras  model  in  the  Beggar  code  used  in  this  study  should  not  be  used  to 
simulate  unsteady  turbulent  flow.  In  every  case,  the  Spalart-Allnraras  model  damped  out 
all  the  unsteady  effects  in  the  flow.  Since  the  DES  model  in  the  Beggar  code  used  Spalart- 
Alllnraras  as  its  subgrid  model,  there  is  concern  that  DES  may  damp  out  the  steady  effects 
as  well.  This  appears  to  have  been  the  case  for  the  oscillating  airfoil  at  aQ  =  11.0°.  The  grid 
may  not  have  been  resolved  fine  enough  in  that  case  though,  as  unsteady  effects  were  not 
damped  out  by  the  DES  model  in  the  turbulent  cavity  simulation  with  the  medium  grid. 
Strong  evidence  was  presented  that  the  Spalart-Allnraras  model  in  Beggar  may  have  been 
written  incorrectly. 

No  conclusions  can  be  made  about  which  turbulence  model  uses  the  least  cpu-time. 
The  largest  effect  that  the  turbulence  model  had  on  average  cpu-time  was  determining  if 
there  should  be  vorticies  or  not.  The  largest  piece  of  cpu-time  is  dedicated  to  converging 
the  inner  iterations,  and  unsteady  phenomena  such  as  vorticies  increase  the  difficulty  of 
converging  these  inner  iterations. 
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V.  GENERIC  STORE  SIMULATION 


5. 1  Case  Presentation 

Chapters  III  and  IV  provide  much  insight  into  the  Beggar  code  and  how  the  turbulence 
models  compare  with  experimental  and  theoretical  data,  but  the  cases  themselves  are  not 
representative  of  the  simulations  performed  day  to  day  by  the  CAT.  To  see  how  the  models 
react  to  a  more  applicable  simulation,  a  MK-84  store  jettison  was  simulated. 

The  MK-84  is  initially  in  an  X-configuration  at  a0  =  —10.0°,  as  shown  in  Figure  66. 
The  MK-84  block  is  broken  into  eight  equal  circumferentially  divided  grids  of  dimension 
138x26x52.  The  store  is  160  inches  long  with  a  maximum  width  of  24.5  inches.  The  initial 
spacing  off  the  store  surface  is  3xl0~5  inches.  The  store  block  resides  inside  an  overlap  grid 
of  dimensions  70x51x45  (750x600x600  inches)  which  is  fixed  to  the  store.  The  store  block 
and  overlap  block  are  both  released  in  the  center  of  a  stationary  cartesian  grid  of  dimensions 
44x43x43  (2200x2000x2000  inches). 

Initial  conditions  are  shown  in  Table  13.  To  achieve  time  accurate  results,  three 
Newton  dt-iterations  were  run  every  time  step.  Each  Newton  dt-iteration  will  solve  80 
symmetric  Gauss-Seidel  inner  iterations  unless  an  inner  convergence  tolerance  of  1  •  10  is 
met. 


Table  13:  Store  jettison  initial  conditions 


Property 

Value 

M 

0.95 

Re 

4.9  •  10&4. 

£qo 

1097.6^ 

Solver  options: 

Steger- Warming  Jacobians 

2nd  Order  Spatial  Differencing 
Roe  Right  Hand  Side 

3  dt-iterations 

The  store  jettison  was  simulated  with  the  Baldwin-Lomax,  Spalart-Allmaras,  and  DES 
models.  Recording  the  change  in  orientation  angle  of  the  store,  the  effect  the  turbulence 
model  has  on  the  store  can  be  recorded.  Turbulent  effects  will  be  greatest  when  the  store  has 
little  inertia.  For  this  reason,  the  simulation  was  initialized  with  no  inertia  in  the  y-direction. 
A  small  initial  time  step  of  At  =  7.6xl0~7  seconds  was  used  to  allow  the  turbulence  model 
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Figure  66:  Initial  Store  Orientation 

to  properly  influence  the  store  trajectory  before  the  inertia  starts  to  dominate.  The  time 
step  used  with  respect  to  iteration  is  shown  in  Table  14.  After  14,000  iterations,  the  Spalart- 
Allmaras  and  DES  models  predicted  negative  densities  so  the  simulation  was  stopped  there. 
It  should  be  noted  that  Baldwin-Lonrax  would  have  continued. 

Table  14:  Store  jettison  time  step  per  iteration 


Iteration  Range 

Time  Step  (sec) 

1-10000 

7.6xlO Y 

10001-11000 

7.6xl0 ti 

11001-13000 

7.6xl0-5 

13001-14000 

1.14xl0~4 

5.2  Results 

Figures  67  to  69  show  the  time  history  of  roll,  pitch,  and  yaw  respectively.  All  of  the 
models  are  in  good  agreement  for  pitch  and  yaw.  Divergence  is  present  in  the  roll  time 
history  plot.  Baldwin-Lomax  predicts  negative  roll  angles  while  Spalart-Allmaras  and  DES 
predict  positive  values.  The  separation  between  Spalart-Allmaras  and  DES  roll  predictions 
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also  appears  to  be  diverging.  This  is  probably  due  to  the  turbulence  models  effect  on  the 
fins. 

Figures  70,  71,  and  72  show  total  pressure  on  and  around  the  store  body.  All  three  fig¬ 
ures  show  the  same  distribution  across  the  store  body.  Figure  73  shows  pressure  coefficients 
along  the  store  for  the  three  models.  All  three  models  predict  the  same  pressure  coefficient 
distribution  except  for  minor  differences  near  the  fins,  showing  that  the  turbulence  model 
will  not  have  much  effect  on  the  motion  of  the  store. 

Figure  74  shows  the  pressure  coefficient  along  the  mid-span  of  one  of  the  upper  fins. 
The  Cp  on  the  upper  surface  of  the  fin  is  the  same  for  all  three  models.  DES  predicts  about 
a  20%  increase  in  Cp  on  the  lower  surface  near  the  quarter  chord.  This  difference  will  not 
create  a  moment  large  enough  to  have  a  substantial  influence  on  a  2000  pound  bomb,  but 
the  fact  that  the  difference  exists  should  be  noted. 

All  runs  for  the  store  jettison  simulation  were  performed  on  5  nodes  of  an  Opteron 
cluster  with  2  processors  per  node  and  four  gigabytes  of  memory  per  node.  The  wall  clock 
times  for  each  run  were  recorded.  The  average  time  per  iteration  for  each  turbulence  model 
is  shown  in  Table  15,  which  shows  the  DES  model  taking  about  |  the  time  of  Spalart- 
Allmaras  model.  This  time  difference  can  be  attributed  to  the  turbulent  PDE  converging 
fast  and  using  less  Gauss-Seidel  inner  iterations  in  the  DES  case.  Baldwin-Lomax  took 
about  1.6  times  longer  to  solve  each  time  step  than  DES. 


Table  15: 


Wa 


1  clock  times  for  various  turbulence  models  for  the  store  jettison  case 


Model 

Time  per  Iteration 
Normalized  by  Baldwin-Lomax 

Baldwin-Lomax 

1.0 

Spalart-Allmaras 

0.81 

DES 

0.61 
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TOT  F 
1.506 


Figure  70:  Pressure  at  t  =  7.62x10  3  sec  -  Baldwin-Lomax 
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0.780 
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Figure  71: 


Pressure  at  t  =  7.62x10  3  sec  -  Spalart-Allmaras 


TOT  F 
1.506 


Figure  72:  Pressure  at  t  =  7.62x10  3  sec  -  DES 


Figure  73:  Pressure  Coefficient  Along  Store  Body 
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VI.  SUMMARY  AND  RECOMMENDATIONS 


6. 1  Summary 

Reviewing  Chapter  III,  conclusions  were  drawn  about  which  solver  options  should  be 
used  while  attempting  to  model  unsteady  turbulent  flow.  To  properly  simulate  unsteady 
turbulent  flow,  it  is  important  to  reduce  artificial  dissipation  and  accurately  calculate  each 
time  step.  The  inviscid  convecting  vortex  provided  much  useful  data  about  how  each  solver 
option  contributed  to  artificial  dissipation,  while  the  shock-tube  and  2D  cylinder  were  most 
beneficial  in  the  area  of  time  step  convergence.  Reviewing  the  non-turbulent  cases,  the 
conclusions  are  drawn  that  the  following  solver  options  should  be  implemented  in  the  Beggar 
code  to  simulate  unsteady  turbulent  flow: 

1.  A  second  order  temporal  discretization 

2.  Between  three  and  five  Newton  dt-iterations 

Considering  all  three  cases  in  the  turbulent  study  in  Chapter  IV,  it  is  determined 
that  the  Spalart-Allmaras  model  as  currently  implemented  should  not  be  used  to  simulate 
unsteady  turbulent  flow.  In  every  case,  the  Spalart-Allmaras  model  damped  out  most  of  the 
unsteady  effects  in  the  flow.  Since  the  DES  model  in  the  Beggar  code  used  Spalart-Allmaras 
as  its  subgrid  model,  there  is  concern  that  DES  may  damp  out  the  unsteady  effects  as  well. 
This  appears  to  have  been  the  case  for  the  oscillating  airfoil  at  aQ  =  11.0°.  However,  the 
grid  may  not  have  been  resolved  fine  enough,  as  unsteady  effects  were  not  damped  out  by 
the  DES  model  in  the  turbulent  cavity  simulation  with  the  medium  grid.  Unfortunately, 
time  and  computational  constraints  prevented  further  study  of  the  DES  model  on  grids 
fine  enough  where  DES  performs  best.  Therefore  a  conclusion  on  whether  the  DES  model 
should  be  used  to  model  unsteady  turbulent  flow  cannot  be  made. 

A  strong  case  has  been  made  that  the  Spalart-Allmaras  model  may  not  have  been 
written  correctly  into  Beggar.  (17:Murman)  notes  that  the  Spalart-Allmaras  model  produces 
excessively  large  values  of  turbulent  eddy- viscosity  in  the  region  of  off-surface  vorticies. 
This  property  of  the  model  may  be  washing  away  the  unsteady  effects,  but  the  fact  that 
the  Beggar  version  of  Spalart-Allmaras  seems  to  be  damping  unsteady  effects  more  than 
other  codes  leads  to  the  conclusion  that  the  Beggar  version  of  Spalart-Allmaras  needs  to 
be  evaluated.  Recalling  that  the  only  difference  between  the  Spalart-Allmaras  and  DES 
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models  is  the  distance  to  the  wall  variable  and  because  DES  provided  much  better  results 
with  the  turbulent  cavity  case,  the  distance  to  the  wall  is  a  likely  suspect  as  to  what  is  being 
improperly  calculated. 

Finally  considering  the  store  jettison  simulation  in  Chapter  V,  the  conclusion  is  drawn 
that  the  choice  of  turbulence  model  does  not  play  a  large  part  in  the  trajectory  of  a  massive 
object  traveling  at  high  velocities.  The  fact  that  the  Baldwin-Lomax  model  took  1.6  times 
longer  to  solve  the  problem  may  be  a  good  reason  to  cease  using  it,  especially  in  steady 
problems  where  Spalart-Allmaras  will  have  the  added  bonus  of  damping  unsteady  effects  in 
the  flow. 

6.2  Recommendations 

This  study  compared  different  turbulence  models  using  the  same  grids.  While  this 
provided  information  about  how  a  turbulence  model  generally  performs,  grid  optimization 
would  provide  additional  information  on  solution  clock  times  and  accuracy.  Each  model 
would  be  optimized  with  differently  dimensioned  grids.  Therefore,  a  good  study  to  run 
would  be  to  optimize  the  grids  for  the  models  being  compared  and  then  compare  clock 
times  and  solutions  to  determine  which  model  best  suits  the  user’s  needs. 

The  instabilities  that  second  order  Roe  combined  with  the  three  point  backward  time 
temporal  discretization  showed  in  the  inviscid  vortex  case  should  be  investigated. 
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Minimum  Pressure 


Appendix  A.  Results  from  Non-Turbulent  Cases 


A.l  Inviscid  Convecting  Vortex 
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Figure  81:  Cycle  vs  Ax  -  Roe  -  Coarse  Grid  -  2nd  Order  Temporal 
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Figure  84:  Cycle  vs  Ay  -  Steger-Warming  -  Coarse  Grid  -  1st  Order  Temporal 
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Figure  88:  Cycle  vs  Ay  -  Roe  -  Fine  Grid  -  1st  Order  Temporal 
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Figure  89:  Cycle  vs  Ax  -  Roe  -  Fine  Grid  -  2nd  Order  Temporal 
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Figure  90:  Cycle  vs  Ay  -  Roe  -  Fine  Grid  -  2nd  Order  Temporal 
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Figure  92:  Cycle  vs  A y  -  Steger- Warming  -  Fine  -  1st  Order  Temporal 
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Figure  95:  Cycle  vs  Ax  -  Roe  -  Coarse  Grid  -  1st  Order  Temporal  -  Halved  Time  Step 
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Figure  96:  Cycle  vs  Ay  -  Roe  -  Coarse  Grid  -  1st  Order  Temporal  -  Halved  Time  Step 
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Figure  98:  Cycle  vs  Ay  -  Roe  -  Coarse  Grid  -  2nd  Order  Temporal  -  Halved  Time  Step 
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Figure  99:  Cycle  vs  Ax  -  Steger- Warming  -  Coarse  Grid  1st  Order  Temporal  -  Halved 

Time  Step 
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Figure  100:  Cycle  vs  Ay  -  Steger- Warming  -  Coarse  Grid  -  1st  Order  Temporal  -  Halved 
Time  Step 
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Figure  101:  Cycle  vs  Ax  -  Steger- Warming  -  Coarse  Grid  -  2nd  Order  Temporal  -  Halved 

Time  Step 


0  1  2  3  4  5  6 

Cycle 

Figure  102:  Cycle  vs  Ay  -  Steger- Warming  -  Coarse  Grid  -  2nd  Order  Temporal  -  Halved 

Time  Step 
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Figure  106:  Cycle  vs  Ay  -  Roe  -  Fine  Grid  -  2nd  Order  Temporal  -  Halved  Time  Step 
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Figure  107:  Cycle  vs  Ax  -  Steger-Warming  -  Fine  Grid  -  1st  Order  Temporal  -  Halved 
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Figure  108: 


Cycle  vs  Ay  -  Steger-Warming  -  Fine  -  1st  Order  Temporal  -  Halved  Time 
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Figure  110:  Cycle  vs  Ay  -  Steger-Warming  -  Fine  Grid  -  2nd  Order  Temporal  -  Halved 

Time  Step 


2D  Vortex  Shedding  Cylinder 


A. 2 


Figure  111:  PSD  of  St  for  Lift  -  Roe  -  1st  Order  Temporal 
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Figure  112:  PSD  of  St  for  Lift  -  Roe  -  2nd  Order  Temporal 


Figure  113:  PSD  of  St  for  Lift  -  Steger- Warming  -  1st  Order  Temporal 
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Figure  114:  PSD  of  St  for  Lift  -  Steger- Warming  -  2nd  Order  Temporal 


Figure  115:  Average  Peak  Drag  vs  dt-iteration 
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Figure  116: 


Average  Drag  vs  dt-iteration 
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Figure  117: 


Average  Peak  Lift  vs  dt-iteration 
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Figure  118:  Average  Strouhal  Number  from  Lift  vs  dt-iteration 
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Appendix  B.  Results  fro?n  Turbulent  Cases 


B.l  Turbulent  Cavity 
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Figure  119:  Sound  Pressure  Level  Spectrum  at  the  K16  Location  -  Coarse  Grid 
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Figure  120:  Sound  Pressure  Level  Spectrum  at  the  K16  Location  -  Medium  Grid 
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Figure  121:  Sound  Pressure  Level  Spectrum  at  the  K16  Location  -  Fine  Grid 
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Figure  123: 


Sound  Pressure  Level  Spectrum  at  the  K18  Location  -  Coarse  Grid 


Sound  Pressure  Level  Spectrum  at  the  K18  Location  -  Medium  Grid 
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Figure  124:  Sound  Pressure  Level  Spectrum  at  the  K18  Location  -  Fine  Grid 
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Figure  125:  Vorticity  -  Baldwin-Lomax  -  Medium  Cavity  Grid 
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Figure  126:  Vorticity  -  Spalart-Allmaras  -  Medium  Cavity  Grid 
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Figure  127:  Vorticity  -  DES  -  Medium  Cavity  Grid 
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Figure  128:  Vorticity  -  Baldwin-Lomax  -  Fine  Cavity  Grid 
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Figure  129:  Vorticity  -  Spalart-Allmaras  -  Fine  Cavity  Grid 
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Figure  130: 


Vorticity  -  Baldwin-Lomax  -  Fine  Cavity  Grid  -  Front 
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Figure  131:  Vorticity  -  Spalart-Allmaras  -  Fine  Cavity  Grid  -  Front 
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